Biostatistical Analysis of Ethanol Pharmacokinetics

Jones_1984 <- read.csv("data/raw/Jones-1984.csv")
Ammon_1996 <- read.csv("data/raw/Ammon-1996.csv")
Yelland_2008 <- read.csv("data/raw/Yelland-2008.csv")
Emberger_2023 <- read.csv("data/raw/Thierauf-Emberger-2023.csv")
Wilkinson_IV <- read.csv("data/raw/Wilkinson-1976-IV.csv")
Wilkinson_PO <- read.csv("data/raw/Wilkinson-1976-PO.csv")
Vestal_1977_elderly <- read.csv("data/raw/Vestal-1977-elderly.csv")
Vestal_1977_young <- read.csv("data/raw/Vestal-1977-young.csv")

Kozawa_2007 <- read.csv("data/raw/Kozawa-2007.csv")
Eriksson_1977 <- read.csv("data/raw/Eriksson-1977.csv")

standard_theme <- theme_bw() +
  theme(plot.title = element_text(size = 30, hjust = 0.5),
    plot.subtitle = element_text(size = 25, hjust = 0.5),
    plot.caption = element_text(size = 20),
    axis.title = element_text(size = 25),
    axis.text.y = element_text(size = 20),
    axis.text.x = element_text(size = 20),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 22),
    strip.text = element_text(size = 22))
theme_set(standard_theme)
set_flextable_defaults(table.layout = "autofit")
flex_default_fun <- function(x) {
  x %>% 
  theme_box() %>%
  font(fontname = "Times New Roman", part = "all") %>% 
  fontsize(size = 13, part = "header") %>% 
  fontsize(size = 12, part = "body") %>% 
  flextable::align(
  i = 1,
  j = NULL,
  align = "center",
  part = "header")
}

1. Human data

1.1 Data unification

Ammon_1996_clean <- Ammon_1996 %>% 
  mutate(case = paste(subject_id, occasion, sep = "-")) %>% 
  mutate(route = route %>% str_remove_all(" infusion")) %>% 
  rename(infusion_duration_min = duration_inf_min)
rm(Ammon_1996)

Jones_1984_clean <- Jones_1984 %>% 
  mutate(subject_id = paste("J", subject_id, sep = "-"),
         sex = case_when(sex == "male" ~ "M",
                         sex == "female" ~ "F"))
rm(Jones_1984)

Yelland_2008_clean <- Yelland_2008 %>% 
  mutate(case = paste(subject_id, occasion, sep = "-"))
rm(Yelland_2008)

Emberger_2023_clean <- Emberger_2023 %>% 
  filter(!is.na(conc_g_per_L)) %>% 
  rename(wt_kg = WT_kg,
         sex = Sex,
         age = Age,
         bh_cm = BH_cm,
         bmi = BMI) %>% 
  mutate(subject_id = str_c("TE", "-", subject_id))
rm(Emberger_2023)

Wilkinson_IV_clean <- Wilkinson_IV %>% 
  filter(!is.na(conc_g_per_L)) %>% 
  mutate(route = route %>% str_remove_all(" infusion")) %>% 
  rename(wt_kg = WT_kg,
         sex = Sex)
rm(Wilkinson_IV)

Wilkinson_PO_clean <- Wilkinson_PO %>% 
  filter(!is.na(conc_g_per_L)) %>% 
  mutate(case = paste(subject_id, occasion, sep = "-")) %>% 
  rename(wt_kg = WT_kg,
         sex = Sex,
         bh_cm = BH_cm,
         bsa_m2 = BSA_m2,
         age = Age)  
rm(Wilkinson_PO)

Vestal_1977_e_pk <- Vestal_1977_elderly %>% 
rename_with(tolower) %>% 
  rename(Cmax = cmax_mg_dl) %>% 
  mutate(Cmax = Cmax / 100,
         route = case_when(
           str_detect(str_to_lower(route), "iv") ~ "IV",
           str_detect(str_to_lower(route), "po") ~ "PO"),
         subject_id = str_c("V-E-", subject_id))
rm(Vestal_1977_elderly)

Vestal_1977_y_pk <- Vestal_1977_young %>% 
rename_with(tolower) %>% 
  rename(Cmax = cmax_mg_dl) %>% 
  mutate(Cmax = Cmax / 100,
         route = case_when(
           str_detect(str_to_lower(route), "iv") ~ "IV",
           str_detect(str_to_lower(route), "po") ~ "PO"),
         subject_id = str_c("V-Y-", subject_id))  
rm(Vestal_1977_young)

1.2. Data plots

1.2.1. Individual curves concentration-time

Ammon_1996_clean %>% 
  select(subject_id, conc_g_per_L, time_min, occasion) %>% 
  mutate(occasion = occasion %>% as.factor) %>% 
  ggplot(aes(color = occasion)) +
  geom_line(aes(x = time_min, y = conc_g_per_L, group = occasion), linewidth = 2) +
  facet_wrap(subject_id~.)+
  scale_x_continuous(breaks = seq(0, 420, by = 50))+
  scale_y_continuous(breaks = seq(0, 1, by = 0.25))+
    labs(title = "Concentration-time plot", subtitle = "Individual curves for each subject", x = "Time, min", y = "Concentration, g/L", color = "Визит", caption = "Source: Ammon-1996")+
  theme(legend.position = "top")

ggplot(Jones_1984_clean, aes(x = time_min, y = conc_g_per_L)) +
  geom_line(aes(group = subject_id)) +
  scale_y_continuous(breaks = seq(0, 1.4, 0.1)) +
  scale_x_continuous(breaks = seq(0, 420, 30)) +
  # coord_cartesian(xlim = c(60, 420)) +
  labs(title = "Concentration-time plot",
       subtitle = "Individual curves for each subject",
       caption = "Source: Jones-1984",
       x = "Time, min",
       y = "Concentration, g/L")

Yelland_2008_clean %>% 
  select(subject_id, conc_g_per_L, time_min, occasion) %>% 
  mutate(occasion = occasion %>% as.factor) %>% 
  ggplot(aes(color = occasion)) +
  geom_line(aes(x = time_min, y = conc_g_per_L, group = occasion), linewidth = 2) +
  facet_wrap(subject_id~.)+
  scale_x_continuous(breaks = seq(0, 420, by = 50))+
  scale_y_continuous(breaks = seq(0, 1, by = 0.25))+
    labs(title = "Concentration-time plot", subtitle = "Individual curves for each subject", x = "Time, min", y = "Concentration, g/L", color = "Occasion", caption = "Source: Yelland-2008")+
  theme(legend.position = "top")

Emberger_2023_clean %>% 
  select(subject_id, conc_g_per_L, time_min) %>% 
  mutate(subject_id = subject_id %>% str_remove_all("TE-")) %>%
  ggplot(aes(color = subject_id)) +
  geom_line(aes(x = time_min, y = conc_g_per_L, group = subject_id), linewidth = 2) +
  scale_x_continuous(breaks = seq(0, 420, by = 50)) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.25)) +
  scale_color_discrete(breaks = as.character(1:10)) +
    labs(title = "Concentration-time plot",
         subtitle = "Individual curves for each subject", 
         x = "Time, min", y = "Concentration, g/L",
         color = "Subject", 
         caption = "Source: Thierauf-Emberger-2023") +
  theme(legend.position = "top")

ggplot(Wilkinson_IV_clean, aes(x = time_min, y = conc_g_per_L)) +
  geom_line(aes(group = subject_id)) +
  scale_y_continuous(breaks = seq(0, 1.4, 0.1)) +
  scale_x_continuous(breaks = seq(0, 510, 30)) +
  # coord_cartesian(xlim = c(60, 420)) +
  labs(title = "Concentration-time plot",
       subtitle = "Individual curves for each subject",
       caption = "Source: Wilkinson-1976-IV",
       x = "Time, min",
       y = "Concentration, g/L")

Wilkinson_PO_clean %>%
  select(subject_id, conc_g_per_L, time_min, occasion) %>% 
  mutate(occasion = occasion %>% as.factor) %>% 
  ggplot(aes(color = occasion)) +
  geom_line(aes(x = time_min, y = conc_g_per_L, group = occasion), linewidth = 2) +
  facet_wrap(subject_id~.)+
  scale_x_continuous(breaks = seq(0, 420, by = 50))+
  scale_y_continuous(breaks = seq(0, 1, by = 0.25))+
    labs(title = "Concentration-time plot", subtitle = "Individual curves for each subject", x = "Time, min", y = "Concentration, g/L", color = "Occasion", caption = "Source: Wilkinson-1976-PO")+
  theme(legend.position = "top")

1.2.2. Mean и standard deviation of concentration

In each time point time is averaged.

Ammon_1996_clean %>% 
  group_by(case) %>%
  mutate(point_id = row_number() %>% as.factor) %>% 
  ungroup() %>%
  mutate(occasion = occasion %>% as.factor) %>% 
  group_by(occasion, point_id) %>% 
  summarise(conc_avg = mean(conc_g_per_L), conc_max = mean(conc_g_per_L) + sd(conc_g_per_L), conc_min = mean(conc_g_per_L) - sd(conc_g_per_L), time_min_avg = mean(time_min)) %>% 
  ungroup() %>% 
  ggplot() +
  geom_line(aes(x = time_min_avg, y = conc_avg, group = occasion, color = occasion)) +
  geom_errorbar(aes(x = time_min_avg, ymin = conc_min, ymax = conc_max), width = 6, linewidth = 0.5) +
  scale_x_continuous(breaks = seq(0, 400, by = 50)) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.25)) +    
  labs(title = "Concentration-time plot",
       subtitle = "Mean ± SD",
       x = "Time, min",
       y = "Concentration, g/L",
       color = "Occasion",
       caption = "Source: Ammon-1996")

Jones_1984_clean %>% 
  group_by(subject_id) %>%
  mutate(point_id = row_number() %>% as.factor) %>% 
  ungroup() %>% 
  group_by(point_id) %>% 
  summarise(conc_avg = mean(conc_g_per_L), conc_max = mean(conc_g_per_L) + sd(conc_g_per_L), conc_min = mean(conc_g_per_L) - sd(conc_g_per_L), time_min_avg = mean(time_min)) %>% 
  ggplot() +
  geom_line(aes(x = time_min_avg, y = conc_avg)) +
  geom_errorbar(aes(x = time_min_avg, ymin = conc_min, ymax = conc_max), width = 6, linewidth = 0.5) +
  scale_x_continuous(breaks = seq(0, 400, by = 50)) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.25)) +    
  labs(title = "Concentration-time plot",
       subtitle = "Mean ± SD",
       x = "Time, min",
       y = "Concentration, g/L",
       caption = "Source: Jones-1984")

Yelland_2008_clean %>%
  group_by(case) %>%
  mutate(point_id = row_number() %>% as.factor) %>% 
  ungroup() %>%
  mutate(occasion = occasion %>% as.factor) %>% 
  group_by(occasion, point_id) %>% 
  summarise(conc_avg = mean(conc_g_per_L), conc_max = mean(conc_g_per_L) + sd(conc_g_per_L), conc_min = mean(conc_g_per_L) - sd(conc_g_per_L), time_min_avg = mean(time_min)) %>% 
  ungroup() %>% 
  ggplot() +
  geom_line(aes(x = time_min_avg, y = conc_avg, group = occasion, color = occasion)) +
  geom_errorbar(aes(x = time_min_avg, ymin = conc_min, ymax = conc_max), width = 6, linewidth = 0.5) +
  scale_x_continuous(breaks = seq(0, 400, by = 50)) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.25)) +    
  labs(title = "Concentration-time plot",
       subtitle = "Mean ± SD",
       x = "Time, min",
       y = "Concentration, g/L",
       color = "Occasion",
       caption = "Source: Yelland-2008")

Emberger_2023_clean %>% 
  group_by(subject_id) %>%
  mutate(point_id = row_number() %>% as.factor) %>% 
  ungroup() %>% 
  group_by(point_id) %>% 
  summarise(conc_avg = mean(conc_g_per_L), conc_max = mean(conc_g_per_L) + sd(conc_g_per_L), conc_min = mean(conc_g_per_L) - sd(conc_g_per_L), time_min_avg = mean(time_min)) %>% 
  ggplot() +
  geom_line(aes(x = time_min_avg, y = conc_avg)) +
  geom_errorbar(aes(x = time_min_avg, ymin = conc_min, ymax = conc_max), width = 6, linewidth = 0.5) +
  scale_x_continuous(breaks = seq(0, 400, by = 50)) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.25)) +    
  labs(title = "Concentration-time plot",
       subtitle = "Mean ± SD",
       x = "Time, min",
       y = "Concentration, g/L",
       caption = "Source: Thierauf-Emberger-2023")

Wilkinson_IV_clean %>% 
  group_by(subject_id) %>%
  mutate(point_id = row_number() %>% as.factor) %>% 
  ungroup() %>% 
  group_by(point_id) %>% 
  summarise(conc_avg = mean(conc_g_per_L), conc_max = mean(conc_g_per_L) + sd(conc_g_per_L), conc_min = mean(conc_g_per_L) - sd(conc_g_per_L), time_min_avg = mean(time_min)) %>% 
  filter(as.numeric(point_id) < 27) %>% 
  ggplot() +
  geom_line(aes(x = time_min_avg, y = conc_avg)) +
  geom_errorbar(aes(x = time_min_avg, ymin = conc_min, ymax = conc_max), width = 6, linewidth = 0.5) +
  scale_x_continuous(breaks = seq(0, 400, by = 50)) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.25)) +    
  labs(title = "Concentration-time plot",
       subtitle = "Mean ± SD",
       x = "Time, min",
       y = "Concentration, g/L",
       caption = "Source: Wilkinson-1976-IV")  

Wilkinson_PO_clean %>%
  group_by(case) %>%
  mutate(point_id = row_number() %>% as.factor) %>% 
  ungroup() %>% 
  group_by(point_id) %>% 
  summarise(conc_avg = mean(conc_g_per_L), conc_max = mean(conc_g_per_L) + sd(conc_g_per_L), conc_min = mean(conc_g_per_L) - sd(conc_g_per_L), time_min_avg = mean(time_min)) %>% 
  ggplot() +
  geom_line(aes(x = time_min_avg, y = conc_avg)) +
  geom_errorbar(aes(x = time_min_avg, ymin = conc_min, ymax = conc_max), width = 6, linewidth = 0.5) +
  scale_x_continuous(breaks = seq(0, 400, by = 50)) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.25)) +    
  labs(title = "Concentration-time plot",
       subtitle = "Mean ± SD",
       x = "Time, min",
       y = "Concentration, g/L",
       caption = "Source: Wilkinson-1976-PO")  

1.3. Evaluation of PK-parameters

#k is how much we add to the tmax index
pharma_fun <- function(conc, t, dose, k = 1) {
  # Modeling for calculation beta, C0, Vd 0 order
  t_sample <- {{t}}[(first(which.max({{conc}})) + k):length({{t}})]
  conc_sample <- conc[(first(which.max({{conc}})) + k):length({{conc}})]
  dose <- first({{dose}})
  model <- lm(conc_sample ~ t_sample)
  Beta <- -as.numeric(model$coefficients[2])
  C_0_0_order <- as.numeric(model$coefficients[1])
  V_d_0_order <- dose / C_0_0_order
  
  # Modeling for calculation kel, C0, Vd 1 order
  # indexes <- conc_sample > 0
  # model_1 <- lm(log(conc_sample[indexes]) ~ t_sample[indexes])
  # Kel <- -as.numeric(model_1$coefficients[2])
  # C_0_1_order <- exp(as.numeric(model_1$coefficients[1]))
  # V_d_1_order <- dose / C_0_1_order
  
  # Calculation of AUC0-t
  
  # if (length({{t}}) < 2) AUC <- NA
  # if (anyNA({{t}}) | anyNA({{conc}})) AUC <- NaN
  t0 = {{t}}[-length({{t}})]
  t1 = {{t}}[-1]
  c0 = {{conc}}[-length({{conc}})]
  c1 = {{conc}}[-1]
  AUC = sum(((c0 + c1)/2)*(t1 - t0))
  
  # Calculation of Cmax
  cmax <- max(conc)
  
  # Calculation of Tmax
  tmax <- t[first(which.max(conc))]
  
  # Calculation of CL (roughly)
  cl = dose/AUC #l/(kg*min)
  # V = CL/Kel 
  
  return(list(Beta = Beta,
              # Kel = Kel,
              C_0_0_order = C_0_0_order,
              V_d_0_order = V_d_0_order,
              # C_0_1_order = C_0_1_order,
              # V_d_1_order = V_d_1_order,
              AUC = AUC,
              Cmax = cmax,
              Tmax = tmax,
              CL = cl))
}

stat_ph <- list(
        "Mean" = ~mean(., na.rm = TRUE) %>% round(4),
        "Median" = ~median(., na.rm = TRUE)%>% round(4),
        "Standard deviation"   = ~sd(., na.rm = TRUE) %>% round(4),
        "Coefficient of variation (CV), %" = ~round((sd(., na.rm = TRUE))/(mean(., na.rm = TRUE)) * 100, 2)
        )

Jones_1984_pk <- Jones_1984_clean %>% 
  arrange(time_min) %>% 
  nest(.by = subject_id, .key = "subject_data") %>% 
  mutate(pharma = purrr::map(subject_data,
                            ~pharma_fun(.$conc_g_per_L,
                                        .$time_min,
                                        .$dose_g_per_kg))) %>% 
  unnest_wider(col = pharma) %>% 
  select(-subject_data)

Jones_1984_pk <- Jones_1984_clean %>% 
  select(-c(time_min, conc_g_per_L)) %>% 
  group_by(subject_id) %>% 
  slice(1) %>% 
  ungroup() %>% 
  right_join(Jones_1984_pk)

## Joining with `by = join_by(subject_id)`

rm(Jones_1984_clean)

Ammon_1996_pk <- Ammon_1996_clean %>% 
  arrange(time_min) %>% 
  nest(.by = case, .key = "subject_data") %>% 
  mutate(pharma = purrr::map(subject_data,
                            ~pharma_fun(.$conc_g_per_L,
                                        .$time_min,
                                        .$dose_g_per_kg))) %>% 
  unnest_wider(col = pharma) %>% 
  select(-subject_data)

Ammon_1996_pk <- Ammon_1996_clean %>%
  select(-c(time_min, conc_g_per_L)) %>%
  group_by(case) %>%
  slice(1) %>%
  ungroup() %>%
  right_join(Ammon_1996_pk, by = join_by(case)) %>% 
  mutate(subject_id = str_c("A-IV-", subject_id)) %>% 
  select(-case)
rm(Ammon_1996_clean)

Yelland_2008_pk <- Yelland_2008_clean %>% 
  arrange(time_min) %>% 
  nest(.by = case, .key = "subject_data") %>% 
  mutate(pharma = purrr::map(subject_data,
                            ~pharma_fun(.$conc_g_per_L,
                                        .$time_min,
                                        .$dose_g_per_kg))) %>% 
  unnest_wider(col = pharma) %>%
  select(-subject_data)

Yelland_2008_pk <- Yelland_2008_clean %>%
  select(-c(time_min, conc_g_per_L)) %>%
  group_by(case) %>%
  slice(1) %>%
  ungroup() %>%
  right_join(Yelland_2008_pk, by = join_by(case)) %>%
  mutate(subject_id = str_c("Y-", subject_id)) %>% 
  select(-case)
rm(Yelland_2008_clean)

Emberger_2023_pk <- Emberger_2023_clean %>% 
  arrange(time_min) %>% 
  nest(.by = subject_id, .key = "subject_data") %>% 
  mutate(pharma = purrr::map(subject_data,
                            ~pharma_fun(.$conc_g_per_L,
                                        .$time_min,
                                        .$dose_g_per_kg))) %>% 
  unnest_wider(col = pharma) %>% 
  select(-subject_data)

Emberger_2023_pk <- Emberger_2023_clean %>% 
  select(-c(time_min, conc_g_per_L)) %>% 
  group_by(subject_id) %>% 
  slice(1) %>% 
  ungroup() %>% 
  right_join(Emberger_2023_pk)

## Joining with `by = join_by(subject_id)`

rm(Emberger_2023_clean)

Wilkinson_IV_pk <- Wilkinson_IV_clean %>% 
  arrange(time_min) %>% 
  nest(.by = subject_id, .key = "subject_data") %>% 
  mutate(pharma = purrr::map(subject_data,
                            ~pharma_fun(.$conc_g_per_L,
                                        .$time_min,
                                        .$dose_g_per_kg))) %>% 
  unnest_wider(col = pharma) %>% 
  select(-subject_data)

Wilkinson_IV_pk <- Wilkinson_IV_clean %>% 
  select(-c(time_min, conc_g_per_L)) %>% 
  group_by(subject_id) %>% 
  slice(1) %>% 
  ungroup() %>% 
  right_join(Wilkinson_IV_pk) %>% 
  mutate(subject_id = str_c("W-IV-", subject_id))

## Joining with `by = join_by(subject_id)`

rm(Wilkinson_IV_clean)

Wilkinson_PO_pk <- Wilkinson_PO_clean %>% 
  arrange(time_min) %>% 
  nest(.by = case, .key = "subject_data") %>% 
  mutate(pharma = purrr::map(subject_data,
                            ~pharma_fun(.$conc_g_per_L,
                                        .$time_min,
                                        .$dose_g_per_kg))) %>% 
  unnest_wider(col = pharma) %>% 
  select(-subject_data)

Wilkinson_PO_pk <- Wilkinson_PO_clean %>%
  select(-c(time_min, conc_g_per_L)) %>%
  group_by(case) %>%
  slice(1) %>%
  ungroup() %>%
  right_join(Wilkinson_PO_pk, by = join_by(case)) %>% 
  mutate(subject_id = str_c("W-PO-", subject_id)) %>% 
  select(-case)
rm(Wilkinson_PO_clean)

1.4. Binding datasets

main_df <- bind_rows(Ammon_1996_pk,
                     Emberger_2023_pk,
                     Jones_1984_pk,
                     Wilkinson_IV_pk,
                     Wilkinson_PO_pk,
                     Yelland_2008_pk,
                     Vestal_1977_e_pk,
                     Vestal_1977_y_pk)

main_df <- main_df %>% 
  mutate(across(c(where(is.character), occasion, infusion_duration_min, -subject_id), ~as.factor(.x)))

# Only mean values by visits (the assumption of independence of observations is fulfilled)
WILK <- main_df %>% 
  filter(study_id == "Wilkinson-1976-PO") %>% 
  arrange(desc(dose_g_per_kg), .by_group = TRUE) %>% 
  group_by(subject_id) %>% 
  slice(1)
order_0_df <- main_df %>%
  filter(study_id != "Wilkinson-1976-PO") %>%
  group_by(subject_id) %>% 
  mutate(across(c(Beta, C_0_0_order, V_d_0_order, AUC, Cmax, Tmax), ~mean(.x))) %>%
  slice(1) %>% 
  ungroup() %>% 
  select(-occasion) %>% 
  bind_rows(WILK) %>% 
  select(-c(
    # Kel, 
    occasion #,
    # contains("1_order")
    )) 

# nrow(order_0_df)
# length(unique(order_0_df$subject_id))
# order_1_df <- main_df %>% filter(Cmax < 0.23)

1.5. Exploratory analysis

1.5.1. Tables

main_df %>% 
  select(-subject_id) %>% 
  gtsummary::tbl_summary()

Characteristic	N = 2061
study_id
Ammon_1996	24 (12%)
Jones_1984	48 (23%)
Thierauf_Emberger_2023	10 (4.9%)
Vestal-1977-elderly	25 (12%)
Vestal-1977-young	25 (12%)
Wilkinson-1976-IV	6 (2.9%)
Wilkinson-1976-PO	32 (16%)
Yelland_2008	36 (17%)
occasion
1	32 (35%)
2	32 (35%)
3	20 (22%)
4	8 (8.7%)
Unknown	114
dose_g_per_kg	0.57 (0.52, 0.68)
route
IV	80 (39%)
PO	126 (61%)
infusion_duration_min
60	74 (93%)
120	6 (7.5%)
Unknown	126
food
fasted	136 (66%)
fed	70 (34%)
site
capillary	86 (42%)
venous	120 (58%)
Beta	0.0022 (0.0019, 0.0026)
Unknown	50
C_0_0_order	0.93 (0.60, 1.00)
Unknown	50
V_d_0_order	0.68 (0.60, 0.72)
Unknown	50
AUC	142 (87, 196)
Unknown	50
Cmax	0.85 (0.65, 1.20)
Tmax	60 (43, 75)
Unknown	50
CL	0.0039 (0.0033, 0.0051)
Unknown	50
wt_kg	75 (71, 82)
Unknown	108
sex
F	5 (3.4%)
M	141 (97%)
Unknown	60
age	32 (23, 61)
Unknown	114
bh_cm	178 (172, 179)
Unknown	164
bmi	23.6 (22.3, 25.2)
Unknown	196
bsa_m2	1.89 (1.83, 2.00)
Unknown	124
lbm_kg	56 (51, 61)
Unknown	156
1 n (%); Median (Q1, Q3)

1.  PK-parameters (endpoints)
•   AUC, Cmax, 
•   C0_0, C0_1
•   Beta, Kel, 
•   Tmax
•   Vd_0,Vd_1
•   CL
2. Experimental conditions
•   Dose (g/kg)
•   Route (IV/PO)
•   Occasion (1/2/NA)
•   Food (fed/fasted)
•   Infusion_duration_min
•   Study_id
3.  Patient's personal data
•   WT, BH, BMI, BSA, LBM
•   Sex, Age

label(main_df$study_id)   <- "Study"
label(main_df$sex)   <- "Sex"
label(main_df$dose_g_per_kg)    <- "Dose"
label(main_df$food) <- "Food"
label(main_df$site) <- "Blood sampling"
label(main_df$age) <- "Age"
label(main_df$route) <- "Route"

units(main_df$dose_g_per_kg)   <- "g per kg"
units(main_df$age)   <- "years"

tbl_output <- table1(~ study_id + dose_g_per_kg + food + site + age + sex|route, data = main_df, render.missing = NULL, overall=c(left="Total"))


tbl_output

	Total (N=206)	IV (N=80)	PO (N=126)
Study
Ammon_1996	24 (11.7%)	24 (30.0%)	0 (0%)
Jones_1984	48 (23.3%)	0 (0%)	48 (38.1%)
Thierauf_Emberger_2023	10 (4.9%)	0 (0%)	10 (7.9%)
Vestal-1977-elderly	25 (12.1%)	25 (31.3%)	0 (0%)
Vestal-1977-young	25 (12.1%)	25 (31.3%)	0 (0%)
Wilkinson-1976-IV	6 (2.9%)	6 (7.5%)	0 (0%)
Wilkinson-1976-PO	32 (15.5%)	0 (0%)	32 (25.4%)
Yelland_2008	36 (17.5%)	0 (0%)	36 (28.6%)
Dose (g per kg)
Mean (SD)	0.557 (0.157)	0.490 (0.125)	0.600 (0.160)
Median [Min, Max]	0.570 [0.127, 0.700]	0.570 [0.300, 0.620]	0.680 [0.127, 0.700]
Food
fasted	136 (66.0%)	56 (70.0%)	80 (63.5%)
fed	70 (34.0%)	24 (30.0%)	46 (36.5%)
Blood sampling
capillary	86 (41.7%)	6 (7.5%)	80 (63.5%)
venous	120 (58.3%)	74 (92.5%)	46 (36.5%)
Age (years)
Mean (SD)	42.1 (19.9)	53.3 (17.6)	28.7 (13.0)
Median [Min, Max]	32.3 [20.7, 80.7]	56.8 [20.7, 80.7]	23.0 [21.0, 67.0]
Sex
F	5 (2.4%)	0 (0%)	5 (4.0%)
M	141 (68.4%)	56 (70.0%)	85 (67.5%)

main_df_long <- main_df %>% 
  select(where(is.numeric)) %>% 
  pivot_longer(
    cols = everything(),
    names_to = "variable",
    values_to = "value"
  )

ggplot(main_df_long, aes(x = value)) +
  geom_histogram(bins = 30) +
  facet_wrap(~ variable, scales = "free") +
  theme_bw() +
  labs(x = NULL, y = NULL)

ggplot(main_df_long, aes(x = value)) +
  geom_density(na.rm = TRUE) +
  facet_wrap(~ variable, scales = "free") +
  theme_bw() +
  labs(x = NULL, y = NULL)

We can take the logarithm of variables with right-sided skewness (V_d_0, AUC, CL)

main_df_log_long <- main_df %>% 
  mutate(across(c(
    # Kel,
    V_d_0_order,
    # V_d_1_order,
    # C_0_1_order,
    AUC,
    CL), function(x) log(x))) %>% 
  select(where(is.numeric)) %>% 
  pivot_longer(
    cols = everything(),
    names_to = "variable",
    values_to = "value"
  )

ggplot(main_df_log_long, aes(x = value)) +
  geom_density(na.rm = TRUE) +
  facet_wrap(~ variable, scales = "free") +
  theme_bw() +
  labs(x = NULL, y = NULL)

main_df_fac_long <- main_df %>% 
  select(where(is.factor)) %>% 
  pivot_longer(
    cols = everything(),
    names_to = "variable",
    values_to = "value"
  )

ggplot(main_df_fac_long, aes(x = value)) +
  geom_bar() +
  facet_wrap(~ variable, scales = "free_x") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = NULL, y = NULL)

1.5.2. Boxplots

theme_custom <- theme(
    axis.text.x = element_text(size = 12),
    axis.text.y = element_text(size = 12),
    axis.title = element_text(size = 14),
    plot.title = element_text(size = 20, hjust = 0.5),
    plot.subtitle = element_text(size = 16, hjust = 0.5),
    legend.title = element_text(size = 16),
    legend.text = element_text(size = 14),
    panel.background = element_rect(fill = "white"), 
    panel.grid = element_blank()
  )

gender_colours <- c(
      "M" = "#8694AC",
      "F" = "#f28482"
    )

route_colours <- c(
      "IV" = "#DDB892",
      "PO" = "#BAC893"
    )

food_colours <- c(
      "fed" = "#FFB35C",
      "fasted" = "#CDE8F4"
    )

site_colours <- c(
      "capillary" = "#FE3448",
      "venous" = "#7C7F65"
    )

gender_fill <- scale_fill_manual(
    name = "Sex",
    values = gender_colours
  )

route_fill <- scale_fill_manual(
    name = "Route",
    values = route_colours
  )

food_fill <- scale_fill_manual(
    name = "Food",
    values = food_colours
  )

site_fill <- scale_fill_manual(
    name = "Site",
    values = site_colours
  )


gender_color <- scale_colour_manual(
    name = "Sex",
    values = gender_colours
  )

p1 <- main_df %>% 
  filter(!is.na(sex)) %>% 
  ggplot(aes(x = sex, y = Tmax)) +
  geom_violin(aes(fill = sex)) +
  geom_boxplot(width = 0.1, color = "black", alpha = 0.5)+
  geom_jitter(color = "gray40", alpha = 0.5, size = 1.5)+
  theme(
    legend.position = "none"
  ) + 
  theme_custom + 
  gender_fill

p2 <- main_df %>% 
  filter(!is.na(sex)) %>%
  ggplot()+
  geom_histogram(aes(x = Tmax, fill = sex), colour = "gray40") + 
  theme_custom + 
  gender_fill

p1 + p2 + plot_annotation(title = "Distribution of Tmax by gender", theme = theme(plot.title = element_text(size = 20)))

p3 <- main_df %>% 
  filter(!is.na(sex)) %>%
  ggplot(aes(x = sex, y = Cmax)) +
  geom_violin(aes(fill = sex)) +
  geom_boxplot(width = 0.1, color = "black", alpha = 0.5)+
  geom_jitter(color = "gray40", alpha = 0.5, size = 1.5)+
  theme(
    legend.position = "none"
  ) + 
  theme_custom + 
  gender_fill

p4 <-  main_df %>% 
  filter(!is.na(sex)) %>%
  ggplot()+
  geom_histogram(aes(x = Cmax, fill = sex), colour = "gray40") + 
  theme_custom + 
  gender_fill

p3 + p4 + plot_annotation(title = "Distribution of Cmax by gender", theme = theme(plot.title = element_text(size = 20)))

p5 <- main_df %>% 
  filter(!is.na(sex)) %>% 
  ggplot(aes(x = sex, y = AUC)) +
  geom_violin(aes(fill = sex)) +
  geom_boxplot(width = 0.1, color = "black", alpha = 0.5)+
  geom_jitter(color = "gray40", alpha = 0.5, size = 1.5)+
  theme(
    legend.position = "none"
  ) + 
  theme_custom + 
  gender_fill

p6 <- main_df %>% 
  filter(!is.na(sex)) %>%
  ggplot()+
  geom_histogram(aes(x = AUC, fill = sex), colour = "gray40") + 
  theme_custom + 
  gender_fill

p5 + p6 + plot_annotation(title = "Distribution of AUC by gender", theme = theme(plot.title = element_text(size = 20)))

p7 <- main_df %>% 
  ggplot(aes(x = route, y = Tmax)) +
  geom_violin(aes(fill = route)) +
  geom_boxplot(width = 0.1, color = "black", alpha = 0.5)+
  geom_jitter(color = "gray40", alpha = 0.5, size = 1.5)+
  theme(
    legend.position = "none"
  ) + 
  theme_custom +
  route_fill

p8 <- main_df %>% 
  ggplot()+
  geom_histogram(aes(x = Tmax, fill = route), colour = "gray40") + 
  theme_custom +
  route_fill

p7 + p8 + plot_annotation(title = "Distribution of Tmax by route", theme = theme(plot.title = element_text(size = 20)))

p9 <- main_df %>% 
  ggplot(aes(x = route, y = Cmax)) +
  geom_violin(aes(fill = route)) +
  geom_boxplot(width = 0.1, color = "black", alpha = 0.5)+
  geom_jitter(color = "gray40", alpha = 0.5, size = 1.5)+
  theme(
    legend.position = "none"
  ) + 
  theme_custom + 
  route_fill

p10 <-  main_df %>% 
  ggplot()+
  geom_histogram(aes(x = Cmax, fill = route), colour = "gray40") + 
  theme_custom + 
  route_fill

p9 + p10 + plot_annotation(title = "Distribution of Cmax by route", theme = theme(plot.title = element_text(size = 20)))

p11 <- main_df %>% 
  ggplot(aes(x = route, y = AUC)) +
  geom_violin(aes(fill = route)) +
  geom_boxplot(width = 0.1, color = "black", alpha = 0.5)+
  geom_jitter(color = "gray40", alpha = 0.5, size = 1.5)+
  theme(
    legend.position = "none"
  ) + 
  theme_custom + 
  route_fill

p12 <- main_df %>% 
  filter(!is.na(sex)) %>%
  ggplot()+
  geom_histogram(aes(x = AUC, fill = route), colour = "gray40") + 
  theme_custom + 
  route_fill

p11 + p12 + plot_annotation(title = "Distribution of AUC by route", theme = theme(plot.title = element_text(size = 20)))

p13 <- main_df %>% 
  ggplot(aes(x = food, y = Tmax)) +
  geom_violin(aes(fill = food)) +
  geom_boxplot(width = 0.1, color = "black", alpha = 0.5)+
  geom_jitter(color = "gray40", alpha = 0.5, size = 1.5)+
  theme(
    legend.position = "none"
  ) + 
  theme_custom +
  food_fill

p14 <- main_df %>% 
  ggplot()+
  geom_histogram(aes(x = Tmax, fill = food), colour = "gray40") + 
  theme_custom +
  food_fill

p13 + p14 + plot_annotation(title = "Distribution of Tmax by food consumption", theme = theme(plot.title = element_text(size = 20)))

p15 <- main_df %>% 
  ggplot(aes(x = food, y = Cmax)) +
  geom_violin(aes(fill = food)) +
  geom_boxplot(width = 0.1, color = "black", alpha = 0.5)+
  geom_jitter(color = "gray40", alpha = 0.5, size = 1.5)+
  theme(
    legend.position = "none"
  ) + 
  theme_custom + 
  food_fill

p16 <-  main_df %>% 
  ggplot()+
  geom_histogram(aes(x = Cmax, fill = food), colour = "gray40") + 
  theme_custom + 
  food_fill

p15 + p16 + plot_annotation(title = "Distribution of Cmax by food consumption", theme = theme(plot.title = element_text(size = 20)))

p17 <- main_df %>% 
  ggplot(aes(x = food, y = AUC)) +
  geom_violin(aes(fill = food)) +
  geom_boxplot(width = 0.1, color = "black", alpha = 0.5)+
  geom_jitter(color = "gray40", alpha = 0.5, size = 1.5)+
  theme(
    legend.position = "none"
  ) + 
  theme_custom + 
  food_fill

p18 <- main_df %>% 
  filter(!is.na(sex)) %>%
  ggplot()+
  geom_histogram(aes(x = AUC, fill = food), colour = "gray40") + 
  theme_custom + 
  food_fill

p17 + p18 + plot_annotation(title = "Distribution of AUC by food consumption", theme = theme(plot.title = element_text(size = 20)))

p19 <- main_df %>% 
  ggplot(aes(x = site, y = Tmax)) +
  geom_violin(aes(fill = site)) +
  geom_boxplot(width = 0.1, color = "black", alpha = 0.5)+
  geom_jitter(color = "gray40", alpha = 0.5, size = 1.5)+
  theme(
    legend.position = "none"
  ) + 
  theme_custom +
  site_fill

p20 <- main_df %>% 
  ggplot()+
  geom_histogram(aes(x = Tmax, fill = site), colour = "gray40") + 
  theme_custom +
  site_fill

p19 + p20 + plot_annotation(title = "Distribution of Tmax by blood sampling", theme = theme(plot.title = element_text(size = 20)))

p21 <- main_df %>% 
  ggplot(aes(x = site, y = Cmax)) +
  geom_violin(aes(fill = site)) +
  geom_boxplot(width = 0.1, color = "black", alpha = 0.5)+
  geom_jitter(color = "gray40", alpha = 0.5, size = 1.5)+
  theme(
    legend.position = "none"
  ) + 
  theme_custom + 
  site_fill

p22 <-  main_df %>% 
  ggplot()+
  geom_histogram(aes(x = Cmax, fill = site), colour = "gray40") + 
  theme_custom + 
  site_fill

p21 + p22 + plot_annotation(title = "Distribution of Cmax by blood sampling", theme = theme(plot.title = element_text(size = 20)))

p23 <- main_df %>% 
  ggplot(aes(x = site, y = AUC)) +
  geom_violin(aes(fill = site)) +
  geom_boxplot(width = 0.1, color = "black", alpha = 0.5)+
  geom_jitter(color = "gray40", alpha = 0.5, size = 1.5)+
  theme(
    legend.position = "none"
  ) + 
  theme_custom + 
  site_fill

p24 <-  main_df %>% 
  ggplot()+
  geom_histogram(aes(x = AUC, fill = site), colour = "gray40") + 
  theme_custom + 
  site_fill

p23 + p24 + plot_annotation(title = "Distribution of AUC by blood sampling", theme = theme(plot.title = element_text(size = 20)))

pk_num <- main_df %>% 
  select(where(is.numeric)) %>% 
  select(-lbm_kg, -bsa_m2, -bh_cm, -bmi) #лишние убрал


cor_mat <- cor(pk_num, use = "pairwise.complete.obs")

cor_mat %>% 
  ggcorr(., label = TRUE)

Testing hypotheses about differences between groups by category “food” and “site”

H0: mean (fasted) = mean(fed) H1: mean (fasted) != mean(fed)

t.test(log(AUC) ~ food, data = order_0_df) 

## 
##  Welch Two Sample t-test
## 
## data:  log(AUC) by food
## t = 9.009, df = 48.616, p-value = 6.134e-12
## alternative hypothesis: true difference in means between group fasted and group fed is not equal to 0
## 95 percent confidence interval:
##  0.4552847 0.7167806
## sample estimates:
## mean in group fasted    mean in group fed 
##             5.264317             4.678284

t.test(log(Cmax) ~ food, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  log(Cmax) by food
## t = 10.867, df = 84.518, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group fasted and group fed is not equal to 0
## 95 percent confidence interval:
##  0.4542268 0.6576812
## sample estimates:
## mean in group fasted    mean in group fed 
##            0.1532168           -0.4027372

t.test(Tmax ~ food, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  Tmax by food
## t = -1.2786, df = 93.904, p-value = 0.2042
## alternative hypothesis: true difference in means between group fasted and group fed is not equal to 0
## 95 percent confidence interval:
##  -16.571704   3.589098
## sample estimates:
## mean in group fasted    mean in group fed 
##             62.70968             69.20098

t.test(Beta ~ food, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  Beta by food
## t = -0.18184, df = 39.468, p-value = 0.8566
## alternative hypothesis: true difference in means between group fasted and group fed is not equal to 0
## 95 percent confidence interval:
##  -0.0002730811  0.0002280155
## sample estimates:
## mean in group fasted    mean in group fed 
##          0.002169623          0.002192156

t.test(C_0_0_order ~ food, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  C_0_0_order by food
## t = 4.225, df = 37.009, p-value = 0.0001495
## alternative hypothesis: true difference in means between group fasted and group fed is not equal to 0
## 95 percent confidence interval:
##  0.09480385 0.26952408
## sample estimates:
## mean in group fasted    mean in group fed 
##            0.9811502            0.7989862

t.test(V_d_0_order ~ food, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  V_d_0_order by food
## t = 1.5073, df = 46.389, p-value = 0.1385
## alternative hypothesis: true difference in means between group fasted and group fed is not equal to 0
## 95 percent confidence interval:
##  -0.009793554  0.068237153
## sample estimates:
## mean in group fasted    mean in group fed 
##            0.6772378            0.6480160

H0: mean (capillary) = mean(venous) H1: mean (capillary) != mean(venous)

t.test(log(AUC) ~ site, data = order_0_df) 

## 
##  Welch Two Sample t-test
## 
## data:  log(AUC) by site
## t = 9.009, df = 48.616, p-value = 6.134e-12
## alternative hypothesis: true difference in means between group capillary and group venous is not equal to 0
## 95 percent confidence interval:
##  0.4552847 0.7167806
## sample estimates:
## mean in group capillary    mean in group venous 
##                5.264317                4.678284

t.test(log(Cmax) ~ site, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  log(Cmax) by site
## t = -4.4454, df = 112.65, p-value = 2.065e-05
## alternative hypothesis: true difference in means between group capillary and group venous is not equal to 0
## 95 percent confidence interval:
##  -0.3673048 -0.1408344
## sample estimates:
## mean in group capillary    mean in group venous 
##              -0.1224289               0.1316407

t.test(Tmax ~ site, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  Tmax by site
## t = -1.2786, df = 93.904, p-value = 0.2042
## alternative hypothesis: true difference in means between group capillary and group venous is not equal to 0
## 95 percent confidence interval:
##  -16.571704   3.589098
## sample estimates:
## mean in group capillary    mean in group venous 
##                62.70968                69.20098

t.test(Beta ~ site, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  Beta by site
## t = -0.18184, df = 39.468, p-value = 0.8566
## alternative hypothesis: true difference in means between group capillary and group venous is not equal to 0
## 95 percent confidence interval:
##  -0.0002730811  0.0002280155
## sample estimates:
## mean in group capillary    mean in group venous 
##             0.002169623             0.002192156

t.test(C_0_0_order ~ site, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  C_0_0_order by site
## t = 4.225, df = 37.009, p-value = 0.0001495
## alternative hypothesis: true difference in means between group capillary and group venous is not equal to 0
## 95 percent confidence interval:
##  0.09480385 0.26952408
## sample estimates:
## mean in group capillary    mean in group venous 
##               0.9811502               0.7989862

t.test(V_d_0_order ~ site, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  V_d_0_order by site
## t = 1.5073, df = 46.389, p-value = 0.1385
## alternative hypothesis: true difference in means between group capillary and group venous is not equal to 0
## 95 percent confidence interval:
##  -0.009793554  0.068237153
## sample estimates:
## mean in group capillary    mean in group venous 
##               0.6772378               0.6480160

Testing hypotheses about differences between groups by category “route” and “sex”.

t.test(Beta ~ route, data = order_0_df) 

## 
##  Welch Two Sample t-test
## 
## data:  Beta by route
## t = -3.952, df = 22.053, p-value = 0.0006754
## alternative hypothesis: true difference in means between group IV and group PO is not equal to 0
## 95 percent confidence interval:
##  -0.0007806703 -0.0002433704
## sample estimates:
## mean in group IV mean in group PO 
##      0.001761587      0.002273607

t.test(Beta ~ sex, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  Beta by sex
## t = 4.5937, df = 4.4652, p-value = 0.007747
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  0.0003360822 0.0012656766
## sample estimates:
## mean in group F mean in group M 
##     0.002982604     0.002181725

t.test(V_d_0_order ~ route, data = order_0_df) 

## 
##  Welch Two Sample t-test
## 
## data:  V_d_0_order by route
## t = -4.1909, df = 29.494, p-value = 0.0002314
## alternative hypothesis: true difference in means between group IV and group PO is not equal to 0
## 95 percent confidence interval:
##  -0.10884915 -0.03748701
## sample estimates:
## mean in group IV mean in group PO 
##        0.6074394        0.6806074

t.test(V_d_0_order ~ sex, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  V_d_0_order by sex
## t = -3.2625, df = 4.6626, p-value = 0.02478
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -0.17872778 -0.01926724
## sample estimates:
## mean in group F mean in group M 
##       0.5752923       0.6742898

t.test(log(AUC) ~ route, data = order_0_df) 

## 
##  Welch Two Sample t-test
## 
## data:  log(AUC) by route
## t = -7.8695, df = 25.189, p-value = 3.004e-08
## alternative hypothesis: true difference in means between group IV and group PO is not equal to 0
## 95 percent confidence interval:
##  -0.7934900 -0.4644088
## sample estimates:
## mean in group IV mean in group PO 
##         4.545742         5.174692

t.test(log(AUC) ~ sex, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  log(AUC) by sex
## t = -4.0071, df = 4.3155, p-value = 0.01379
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -1.1037583 -0.2154553
## sample estimates:
## mean in group F mean in group M 
##        4.564367        5.223974

t.test(Cmax ~ route, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  Cmax by route
## t = 7.9694, df = 80.102, p-value = 9.216e-12
## alternative hypothesis: true difference in means between group IV and group PO is not equal to 0
## 95 percent confidence interval:
##  0.3849115 0.6411207
## sample estimates:
## mean in group IV mean in group PO 
##        1.3830588        0.8700427

t.test(Cmax ~ sex, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  Cmax by sex
## t = -5.4861, df = 5.6395, p-value = 0.001876
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -0.7970139 -0.3000288
## sample estimates:
## mean in group F mean in group M 
##        0.668000        1.216521

t.test(Tmax ~ route, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  Tmax by route
## t = 2.5439, df = 24.323, p-value = 0.01772
## alternative hypothesis: true difference in means between group IV and group PO is not equal to 0
## 95 percent confidence interval:
##   3.588895 34.338455
## sample estimates:
## mean in group IV mean in group PO 
##         80.41667         61.45299

t.test(Tmax ~ sex, data = order_0_df)

## 
##  Welch Two Sample t-test
## 
## data:  Tmax by sex
## t = -0.032711, df = 4.5118, p-value = 0.9753
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -43.19236  42.14162
## sample estimates:
## mean in group F mean in group M 
##        63.40000        63.92537

Results: The administration route groups differed in all pharmacokinetic parameters. The gender groups differed in all pharmacokinetic parameters except Tmax.

1.6. Main Regression Analysis

1.6.1. Regression models for AUC

1.6.1.1. Additive model with main variables (dose, route of administration, blood sampling method)

log(AUC) ~ dose_g_per_kg + route + site

fit_AUC_mod <- lm(log(AUC) ~ dose_g_per_kg + route + site, order_0_df) #+food
fit_AUC_submod <- lm(log(AUC) ~ 1, order_0_df)

#broom::tidy(fit_AUC_mod, conf.int = TRUE)
#broom::tidy(fit_AUC_submod, conf.int = TRUE)

summ(fit_AUC_mod, confint = TRUE, digits = 3)

## MODEL INFO:
## Observations: 96 (50 missing obs. deleted)
## Dependent Variable: log(AUC)
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(3,92) = 102.707, p = 0.000
## R² = 0.770
## Adj. R² = 0.763 
## 
## Standard errors:OLS
## ---------------------------------------------------------------
##                         Est.     2.5%    97.5%   t val.       p
## ------------------- -------- -------- -------- -------- -------
## (Intercept)            4.167    3.874    4.461   28.200   0.000
## dose_g_per_kg          1.544    0.934    2.155    5.021   0.000
## routePO                0.085   -0.093    0.262    0.947   0.346
## sitevenous            -0.342   -0.441   -0.243   -6.875   0.000
## ---------------------------------------------------------------

#summary(fit_AUC_mod)
#summary(fit_AUC_submod)

anova(fit_AUC_mod, fit_AUC_submod)

## Analysis of Variance Table
## 
## Model 1: log(AUC) ~ dose_g_per_kg + route + site
## Model 2: log(AUC) ~ 1
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     92  3.3115                                  
## 2     95 14.4023 -3   -11.091 102.71 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_model(fit_AUC_mod)

ggpredict(fit_AUC_mod, terms = c("dose_g_per_kg", "route", "site"), ci_level = 0) %>% plot(show_data = TRUE, log_y = FALSE, jitter = .01) + labs(subtitle = "Adjusted R-squared:  0.7626")

## Model has log-transformed response. Back-transforming predictions to
##   original response scale. Standard errors are still on the transformed
##   scale.

plot_coefs(fit_AUC_mod, omit.coefs = NULL, scale = TRUE, ci_level = 0.95, plot.distributions = FALSE, inner_ci_level = .1, colors = "CUD") + labs(title = "Regression coeffecients with 95% CI")

1.6.1.2. Additive model with additional variables

log(AUC) ~ dose_g_per_kg + route + site + sex

fit_AUC_mod_add <- lm(log(AUC) ~ dose_g_per_kg + route + site + sex, order_0_df) #age cannot be added

summ(fit_AUC_mod, confint = TRUE, digits = 3)

## MODEL INFO:
## Observations: 96 (50 missing obs. deleted)
## Dependent Variable: log(AUC)
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(3,92) = 102.707, p = 0.000
## R² = 0.770
## Adj. R² = 0.763 
## 
## Standard errors:OLS
## ---------------------------------------------------------------
##                         Est.     2.5%    97.5%   t val.       p
## ------------------- -------- -------- -------- -------- -------
## (Intercept)            4.167    3.874    4.461   28.200   0.000
## dose_g_per_kg          1.544    0.934    2.155    5.021   0.000
## routePO                0.085   -0.093    0.262    0.947   0.346
## sitevenous            -0.342   -0.441   -0.243   -6.875   0.000
## ---------------------------------------------------------------

summ(fit_AUC_mod_add, confint = TRUE, digits = 3)

## MODEL INFO:
## Observations: 72 (74 missing obs. deleted)
## Dependent Variable: log(AUC)
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(4,67) = 39.615, p = 0.000
## R² = 0.703
## Adj. R² = 0.685 
## 
## Standard errors:OLS
## ---------------------------------------------------------------
##                         Est.     2.5%    97.5%   t val.       p
## ------------------- -------- -------- -------- -------- -------
## (Intercept)            3.144    2.407    3.880    8.523   0.000
## dose_g_per_kg          3.154    1.797    4.510    4.640   0.000
## routePO                0.116   -0.080    0.311    1.180   0.242
## sitevenous            -0.341   -0.535   -0.147   -3.512   0.001
## sexM                  -0.068   -0.311    0.175   -0.556   0.580
## ---------------------------------------------------------------

# summary(fit_AUC_mod)
# summary(fit_AUC_mod_add)


#anova(fit_AUC_mod, fit_AUC_mod_add) Anova does not work with variable sex

check_model(fit_AUC_mod_add)

ggpredict(fit_AUC_mod_add, terms = c("dose_g_per_kg", "route", "site", "sex")) %>% plot(show_data = TRUE, log_y = F, jitter = .01) #+ labs(subtitle = "Adjusted R-squared:  0.7626")

## Model has log-transformed response. Back-transforming predictions to
##   original response scale. Standard errors are still on the transformed
##   scale.

1.6.1.3. Model with interactions

log(AUC) ~ dose_g_per_kg*route + site

It is logical that AUC will depend on the combination of dose and method of administration.

fit_AUC_mod_int <- lm(log(AUC) ~ dose_g_per_kg*route + site, order_0_df)
# fit_AUC_mod_int2 <- lm(log(AUC) ~ dose_g_per_kg + route*site, order_0_df)
# fit_AUC_mod_int3 <- lm(log(AUC) ~ dose_g_per_kg*site + route, order_0_df)
# fit_AUC_mod_int4 <- lm(log(AUC) ~ dose_g_per_kg*site*route, order_0_df)


summ(fit_AUC_mod_int, confint = TRUE, digits = 3)

## MODEL INFO:
## Observations: 96 (50 missing obs. deleted)
## Dependent Variable: log(AUC)
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(4,91) = 95.676, p = 0.000
## R² = 0.808
## Adj. R² = 0.799 
## 
## Standard errors:OLS
## -----------------------------------------------------------------------
##                                 Est.     2.5%    97.5%   t val.       p
## --------------------------- -------- -------- -------- -------- -------
## (Intercept)                    4.520    4.204    4.837   28.365   0.000
## dose_g_per_kg                  0.672   -0.023    1.367    1.921   0.058
## routePO                       -1.079   -1.649   -0.509   -3.762   0.000
## sitevenous                    -0.358   -0.449   -0.266   -7.791   0.000
## dose_g_per_kg:routePO          2.107    1.118    3.096    4.233   0.000
## -----------------------------------------------------------------------

anova(fit_AUC_mod_int, fit_AUC_submod)

## Analysis of Variance Table
## 
## Model 1: log(AUC) ~ dose_g_per_kg * route + site
## Model 2: log(AUC) ~ 1
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     91  2.7667                                  
## 2     95 14.4023 -4   -11.636 95.676 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_model(fit_AUC_mod_int)

ggpredict(fit_AUC_mod_int, terms = c("dose_g_per_kg", "route", "site"), ci_level = 0) %>% plot(show_data = TRUE, log_y = FALSE, jitter = .01) + labs(subtitle = "Adjusted R-squared:  0.799") + xlab("Dose (g per kg)")

## Model has log-transformed response. Back-transforming predictions to
##   original response scale. Standard errors are still on the transformed
##   scale.

plot_coefs(fit_AUC_mod_int, omit.coefs = NULL, scale = TRUE, ci_level = 0.95, plot.distributions = FALSE, inner_ci_level = .1, colors = "CUD") + labs(title = "Regression coeffecients with 95% CI")

plot_model(fit_AUC_mod_int, show.values = TRUE,  show.intercept = TRUE, value.offset = 0.3,  axis.labels = c("Dose (g per kg)", "Route (PO)", "Blood sampling (venous)", "Dose (g per kg) x Route (PO)", "Intercept"), title = "Regression coeffecients with 95% CI") +
  theme_nice()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the sjPlot package.
##   Please report the issue at <https://github.com/strengejacke/sjPlot/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Coclusions:

Including the interaction between dose_g_per_kg and route increases \(R^2\) by about 3%.

1.6.2. Regression models for Tmax

order_0_df <- order_0_df %>% 
  mutate(bmi = ifelse(is.na(bmi) & !is.na(wt_kg) & !is.na(bh_cm),
                      wt_kg/(bh_cm/100)^2, bmi),
         bsa_m2 = ifelse(is.na(bsa_m2) & !is.na(wt_kg) & !is.na(bh_cm),
                         sqrt((bh_cm * wt_kg)/3600), bsa_m2))

t_max_full_m <- lm(Tmax ~ dose_g_per_kg * food + route + wt_kg, data = order_0_df)

summary(t_max_full_m)

## 
## Call:
## lm(formula = Tmax ~ dose_g_per_kg * food + route + wt_kg, data = order_0_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.900 -10.838   1.341   4.782  42.309 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)   
## (Intercept)            253.9999   102.2243   2.485  0.02303 * 
## dose_g_per_kg         -165.4839   145.4369  -1.138  0.27011   
## foodfed               -247.9988   129.1429  -1.920  0.07080 . 
## routePO                -48.0073    12.3665  -3.882  0.00109 **
## wt_kg                   -0.4782     0.4216  -1.134  0.27153   
## dose_g_per_kg:foodfed  433.3209   221.8195   1.953  0.06648 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.64 on 18 degrees of freedom
##   (122 пропущенных наблюдений удалены)
## Multiple R-squared:  0.6062, Adjusted R-squared:  0.4969 
## F-statistic: 5.542 on 5 and 18 DF,  p-value: 0.002916

check_model(t_max_full_m)

ggResidpanel::resid_xpanel(t_max_full_m)

ggResidpanel::resid_panel(t_max_full_m, plots = c("lev", "cookd"))

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the ggResidpanel package.
##   Please report the issue to the authors.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Estimating model coefficients using sandwich estimator

t_max_coef <- coeftest(t_max_full_m, vcov. = vcovHC, type = "HC3") %>%
broom::tidy(conf.int = TRUE) %>% 
  mutate(term = c("Intercept",
                  "Dose (g/kg)",
                  "Food status: \"Fed\"",
                  "Peroral route",
                  "Weight, kg",
                  "Interaction between\ndose and status \"fed\"")) %>% 
  rename("Variable" = term,
         "β estimate" = estimate,
         "Standard error" = std.error,
         "Wald t-test statistics" = statistic,
         "P-value" = p.value,
         "95%CI low" = conf.low,
         "95%CI high" = conf.high)
t_max_coef %>% 
  mutate(across(where(is.numeric), ~round(.x, 3))) %>% 
  flextable()

Variable	β estimate	Standard error	Wald t-test statistics	P-value	95%CI low	95%CI high
Intercept	254.000	162.510	1.563	0.135	-87.421	595.421
Dose (g/kg)	-165.484	233.769	-0.708	0.488	-656.615	325.647
Food status: "Fed"	-247.999	211.833	-1.171	0.257	-693.044	197.046
Peroral route	-48.007	14.871	-3.228	0.005	-79.249	-16.765
Weight, kg	-0.478	0.675	-0.709	0.488	-1.896	0.939
Interaction between dose and status "fed"	433.321	353.804	1.225	0.236	-309.993	1,176.635

t_max_prediction_plot <- fortify(t_max_full_m)

## Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
## ℹ Please use `broom::augment(<lm>)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ggplot(t_max_prediction_plot, aes(x = .fitted, y = Tmax)) +
  geom_point(size = 2) +
  ylim(c(25,130)) +
  xlim(c(25,130)) +
  geom_abline(slope = 1, intercept = 0, color = "skyblue3",
              linetype = "dashed", linewidth = 2) +
  geom_text(x = 75, y = 130, label = "Adjusted R-squared = 0.497", size = 8) +
  labs(title = "Predicted vs observed Tmax values",
       x = "Fitted values of Tmax, min",
       y = "Observed values of Tmax, min")

t_max_coef %>% 
ggplot(aes(y = `Variable`)) +
  geom_point(aes(x = `β estimate`), size = 3) +
  geom_errorbar(aes(xmin = `95%CI low`, xmax = `95%CI high`), width = 0.5) +
  geom_vline(xintercept = 0, color = "darkred") +
  labs(title = "Coefficients and 95%CI in Tmax model")

1.6.3. Regression models for Cmax

log(Cmax) ~ dose_g_per_kg + route + food + site

order_0_df_factor <- order_0_df %>% 
  mutate(
    route = factor(route),
    food = factor(food),
    site = factor(site)
  )

Cmax_lm2_v1 <- lm(log(Cmax) ~ dose_g_per_kg + route + food + site, data = order_0_df_factor)

summary(Cmax_lm2_v1)

## 
## Call:
## lm(formula = log(Cmax) ~ dose_g_per_kg + route + food + site, 
##     data = order_0_df_factor)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.48366 -0.10779  0.01501  0.10903  0.36264 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.60571    0.12810  -4.728 5.43e-06 ***
## dose_g_per_kg  0.57351    0.26664   2.151   0.0332 *  
## routePO        0.11544    0.07734   1.493   0.1378    
## foodfed       -0.94209    0.07201 -13.083  < 2e-16 ***
## sitevenous     0.77382    0.06083  12.720  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1645 on 141 degrees of freedom
## Multiple R-squared:  0.8369, Adjusted R-squared:  0.8322 
## F-statistic: 180.8 on 4 and 141 DF,  p-value: < 2.2e-16

check_model(Cmax_lm2_v1)

broom::tidy(Cmax_lm2_v1, conf.int = TRUE)

## # A tibble: 5 × 7
##   term          estimate std.error statistic  p.value conf.low conf.high
##   <chr>            <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)     -0.606    0.128      -4.73 5.43e- 6  -0.859     -0.352
## 2 dose_g_per_kg    0.574    0.267       2.15 3.32e- 2   0.0464     1.10 
## 3 routePO          0.115    0.0773      1.49 1.38e- 1  -0.0375     0.268
## 4 foodfed         -0.942    0.0720    -13.1  4.18e-26  -1.08      -0.800
## 5 sitevenous       0.774    0.0608     12.7  3.61e-25   0.654      0.894

# analysis of model linearity using Residuals vs Fitted plot
autoplot(Cmax_lm2_v1, 1, ncol = 1)

# Linearity diagnosis with ggResidpanel library
resid_panel(Cmax_lm2_v1, plots = "resid", smoother = TRUE)

## `geom_smooth()` using formula = 'y ~ x'

# Diagnosis of homoscedasticity of residuals using the Scale-Location plot
autoplot(Cmax_lm2_v1, 3, ncol = 1)

# Diagnosing the normality of residual distribution using a Normal Q-Q plot:
resid_panel(Cmax_lm2_v1, plots = c("qq", "hist"))

# Diagnosis of outliers and anomalous observations using the Residuals vs Leverage plot and Cook's distance plot
resid_panel(Cmax_lm2_v1, plots = c("lev", "cookd"))

log(Cmax) ~ dose_g_per_kg + route + food + site + sex + age

order_0_df_factor <- order_0_df %>% 
  mutate(
    route = factor(route),
    food = factor(food),
    site = factor(site),
    sex = factor(sex)
  )

Cmax_lm2_v2 <- lm(log(Cmax) ~ dose_g_per_kg + route + food + sex + age, data = order_0_df_factor)

summary(Cmax_lm2_v2)

## 
## Call:
## lm(formula = log(Cmax) ~ dose_g_per_kg + route + food + sex + 
##     age, data = order_0_df_factor)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33748 -0.07599 -0.00931  0.11259  0.29027 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.149688   0.490670   0.305  0.76134    
## dose_g_per_kg -0.012206   0.953701  -0.013  0.98983    
## routePO       -0.507512   0.075700  -6.704 6.97e-09 ***
## foodfed       -0.286139   0.093412  -3.063  0.00324 ** 
## sexM           0.087738   0.119204   0.736  0.46449    
## age            0.004959   0.001153   4.301 6.13e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1523 on 62 degrees of freedom
##   (78 пропущенных наблюдений удалены)
## Multiple R-squared:  0.8574, Adjusted R-squared:  0.8459 
## F-statistic: 74.58 on 5 and 62 DF,  p-value: < 2.2e-16

check_model(Cmax_lm2_v1)

broom::tidy(Cmax_lm2_v2, conf.int = TRUE)

## # A tibble: 6 × 7
##   term          estimate std.error statistic       p.value conf.low conf.high
##   <chr>            <dbl>     <dbl>     <dbl>         <dbl>    <dbl>     <dbl>
## 1 (Intercept)    0.150     0.491      0.305  0.761         -0.831     1.13   
## 2 dose_g_per_kg -0.0122    0.954     -0.0128 0.990         -1.92      1.89   
## 3 routePO       -0.508     0.0757    -6.70   0.00000000697 -0.659    -0.356  
## 4 foodfed       -0.286     0.0934    -3.06   0.00324       -0.473    -0.0994 
## 5 sexM           0.0877    0.119      0.736  0.464         -0.151     0.326  
## 6 age            0.00496   0.00115    4.30   0.0000613      0.00265   0.00726

# analysis of model linearity using Residuals vs Fitted plot
autoplot(Cmax_lm2_v2, 1, ncol = 1)

# Linearity diagnosis with ggResidpanel library
resid_panel(Cmax_lm2_v2, plots = "resid", smoother = TRUE)

## `geom_smooth()` using formula = 'y ~ x'

# Diagnosis of homoscedasticity of residuals using the Scale-Location plot
autoplot(Cmax_lm2_v2, 3, ncol = 1)

# Diagnosing the normality of residual distribution using a Normal Q-Q plot:
resid_panel(Cmax_lm2_v2, plots = c("qq", "hist"))

# Diagnosis of outliers and anomalous observations using the Residuals vs Leverage plot and Cook's distance plot
resid_panel(Cmax_lm2_v2, plots = c("lev", "cookd"))

t_max_coef <- coeftest(Cmax_lm2_v2, vcov. = vcovHC, type = "HC3") %>%
broom::tidy(conf.int = TRUE) %>% 
  mutate(term = c("Intercept",
                  "Dose (g/kg)",
                  "Peroral route",
                  "Food status: \"Fed\"",
                  "sex M",
                  "age")) %>% 
  rename("Variable" = term,
         "β estimate" = estimate,
         "Standard error" = std.error,
         "Wald t-test statistics" = statistic,
         "P-value" = p.value,
         "95%CI low" = conf.low,
         "95%CI high" = conf.high)
t_max_coef %>% 
  mutate(across(where(is.numeric), ~round(.x, 3))) %>% 
  flextable()

Variable	β estimate	Standard error	Wald t-test statistics	P-value	95%CI low	95%CI high
Intercept	0.150	1.258	0.119	0.906	-2.365	2.665
Dose (g/kg)	-0.012	2.449	-0.005	0.996	-4.907	4.883
Peroral route	-0.508	0.171	-2.971	0.004	-0.849	-0.166
Food status: "Fed"	-0.286	0.158	-1.813	0.075	-0.602	0.029
sex M	0.088	0.233	0.377	0.708	-0.378	0.553
age	0.005	0.001	4.278	0.000	0.003	0.007

ggpredict(Cmax_lm2_v2, terms = c("age", "route"), ci_level = 0) %>% plot(show_data = TRUE, log_y = FALSE, jitter = .01) + labs(subtitle = "Adjusted R-squared:  0.8459")

## Model has log-transformed response. Back-transforming predictions to
##   original response scale. Standard errors are still on the transformed
##   scale.

plot_model(Cmax_lm2_v2, show.values = TRUE,  show.intercept = TRUE, value.offset = 0.3, digits = 3,title = "Regression coeffecients with 95% CI") +
  geom_hline(yintercept = 0, color = "gray3", size = 0.5) +
  theme_nice()

log(Cmax) ~ age + lbm_kg

Cmax_lm2_v3 <- lm(Cmax ~ age + lbm_kg, data = order_0_df_factor)

summary(Cmax_lm2_v3)

## 
## Call:
## lm(formula = Cmax ~ age + lbm_kg, data = order_0_df_factor)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35699 -0.11888 -0.02071  0.09604  0.42718 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.862906   0.304755   6.113 1.82e-07 ***
## age          0.005261   0.001886   2.789  0.00762 ** 
## lbm_kg      -0.008623   0.004188  -2.059  0.04508 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1935 on 47 degrees of freedom
##   (96 пропущенных наблюдений удалены)
## Multiple R-squared:   0.36,  Adjusted R-squared:  0.3328 
## F-statistic: 13.22 on 2 and 47 DF,  p-value: 2.784e-05

#check_model(Cmax_lm2_v1)

broom::tidy(Cmax_lm2_v3, conf.int = TRUE)

## # A tibble: 3 × 7
##   term        estimate std.error statistic     p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>       <dbl>    <dbl>     <dbl>
## 1 (Intercept)  1.86      0.305        6.11 0.000000182  1.25     2.48    
## 2 age          0.00526   0.00189      2.79 0.00762      0.00147  0.00906 
## 3 lbm_kg      -0.00862   0.00419     -2.06 0.0451      -0.0170  -0.000197

# analysis of model linearity using Residuals vs Fitted plot
autoplot(Cmax_lm2_v3, 1, ncol = 1)

# Linearity diagnosis with ggResidpanel library
resid_panel(Cmax_lm2_v3, plots = "resid", smoother = TRUE)

## `geom_smooth()` using formula = 'y ~ x'

# Diagnosis of homoscedasticity of residuals using the Scale-Location plot
autoplot(Cmax_lm2_v3, 3, ncol = 1)

# Diagnosing the normality of residual distribution using a Normal Q-Q plot:
resid_panel(Cmax_lm2_v3, plots = c("qq", "hist"))

# Diagnosis of outliers and anomalous observations using the Residuals vs Leverage plot and Cook's distance plot
resid_panel(Cmax_lm2_v3, plots = c("lev", "cookd"))

1.6.4. Regression models for C_0

order_0_df

## # A tibble: 146 × 21
##    study_id   subject_id dose_g_per_kg route infusion_duration_min food  site  
##    <fct>      <chr>              <dbl> <fct> <fct>                 <fct> <fct> 
##  1 Ammon_1996 A-IV-1               0.3 IV    60                    fed   venous
##  2 Ammon_1996 A-IV-10              0.3 IV    60                    fed   venous
##  3 Ammon_1996 A-IV-11              0.3 IV    60                    fed   venous
##  4 Ammon_1996 A-IV-12              0.3 IV    60                    fed   venous
##  5 Ammon_1996 A-IV-2               0.3 IV    60                    fed   venous
##  6 Ammon_1996 A-IV-3               0.3 IV    60                    fed   venous
##  7 Ammon_1996 A-IV-4               0.3 IV    60                    fed   venous
##  8 Ammon_1996 A-IV-5               0.3 IV    60                    fed   venous
##  9 Ammon_1996 A-IV-6               0.3 IV    60                    fed   venous
## 10 Ammon_1996 A-IV-7               0.3 IV    60                    fed   venous
## # ℹ 136 more rows
## # ℹ 14 more variables: Beta <dbl>, C_0_0_order <dbl>, V_d_0_order <dbl>,
## #   AUC <dbl>, Cmax <dbl>, Tmax <dbl>, CL <dbl>, wt_kg <dbl>, sex <fct>,
## #   age <dbl>, bh_cm <int>, bmi <dbl>, bsa_m2 <dbl>, lbm_kg <dbl>

C_0_0_lm1_v1 <- lm(log(C_0_0_order) ~ dose_g_per_kg + route + food + site + study_id, data = order_0_df)
summary(C_0_0_lm1_v1)

## 
## Call:
## lm(formula = log(C_0_0_order) ~ dose_g_per_kg + route + food + 
##     site + study_id, data = order_0_df)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.273548 -0.053489 -0.001246  0.066086  0.283716 
## 
## Coefficients: (3 not defined because of singularities)
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    -0.81812    0.11202  -7.303 1.14e-10 ***
## dose_g_per_kg                   0.57296    0.47240   1.213   0.2284    
## routePO                         0.43382    0.19331   2.244   0.0273 *  
## foodfed                        -0.04388    0.06301  -0.696   0.4879    
## sitevenous                           NA         NA      NA       NA    
## study_idJones_1984             -0.02089    0.05113  -0.409   0.6838    
## study_idThierauf_Emberger_2023  0.03276    0.07954   0.412   0.6814    
## study_idWilkinson-1976-IV       0.42749    0.18760   2.279   0.0251 *  
## study_idWilkinson-1976-PO            NA         NA      NA       NA    
## study_idYelland_2008                 NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09982 on 89 degrees of freedom
##   (50 пропущенных наблюдений удалены)
## Multiple R-squared:  0.8394, Adjusted R-squared:  0.8286 
## F-statistic: 77.52 on 6 and 89 DF,  p-value: < 2.2e-16

check_model(C_0_0_lm1_v1)

Adjusted R-squared: 0.8286

1. Posterior Predictive Check

The density lines of the models and observations almost coincide. This means that the model adequately reproduces the distribution of C_0_0.

2. Linearity

The green smoothed line is close to horizontal, but there is a slight systematic trend visible at the right end (with large fitted values). There is no significant violation of linearity. It is possible that there is insufficient interaction or nonlinear effect of one of the predictors.

3. Homogeneity of variance There is no heteroscedasticity.

4. Influential observations

Several influential observations are visible, marked with numbers (78, 92, 62, 95, 67)

5. Normality of residuals

The Q–Q chart looks acceptable.

6. Collinearty VIF dose_g_per_kg, route, study_id - High ≥ 10

I will try to remove unnecessary indicators, for example, in each study, 1 level of food.

C_0_0_lm1_v2 <- lm(log(C_0_0_order) ~ study_id, data = order_0_df)
summary(C_0_0_lm1_v2)

## 
## Call:
## lm(formula = log(C_0_0_order) ~ study_id, data = order_0_df)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.260871 -0.053489  0.001465  0.066086  0.283716 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    -0.69011    0.02889  -23.89   <2e-16 ***
## study_idJones_1984              0.67454    0.03230   20.88   <2e-16 ***
## study_idThierauf_Emberger_2023  0.61441    0.04285   14.34   <2e-16 ***
## study_idWilkinson-1976-IV       0.63085    0.05004   12.61   <2e-16 ***
## study_idWilkinson-1976-PO       0.65411    0.04568   14.32   <2e-16 ***
## study_idYelland_2008            0.66301    0.04086   16.23   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1001 on 90 degrees of freedom
##   (50 пропущенных наблюдений удалены)
## Multiple R-squared:  0.8367, Adjusted R-squared:  0.8277 
## F-statistic: 92.25 on 5 and 90 DF,  p-value: < 2.2e-16

check_model(C_0_0_lm1_v2)

Adjusted R-squared: 0.8277 Removed collinearity.

order_0_df

## # A tibble: 146 × 21
##    study_id   subject_id dose_g_per_kg route infusion_duration_min food  site  
##    <fct>      <chr>              <dbl> <fct> <fct>                 <fct> <fct> 
##  1 Ammon_1996 A-IV-1               0.3 IV    60                    fed   venous
##  2 Ammon_1996 A-IV-10              0.3 IV    60                    fed   venous
##  3 Ammon_1996 A-IV-11              0.3 IV    60                    fed   venous
##  4 Ammon_1996 A-IV-12              0.3 IV    60                    fed   venous
##  5 Ammon_1996 A-IV-2               0.3 IV    60                    fed   venous
##  6 Ammon_1996 A-IV-3               0.3 IV    60                    fed   venous
##  7 Ammon_1996 A-IV-4               0.3 IV    60                    fed   venous
##  8 Ammon_1996 A-IV-5               0.3 IV    60                    fed   venous
##  9 Ammon_1996 A-IV-6               0.3 IV    60                    fed   venous
## 10 Ammon_1996 A-IV-7               0.3 IV    60                    fed   venous
## # ℹ 136 more rows
## # ℹ 14 more variables: Beta <dbl>, C_0_0_order <dbl>, V_d_0_order <dbl>,
## #   AUC <dbl>, Cmax <dbl>, Tmax <dbl>, CL <dbl>, wt_kg <dbl>, sex <fct>,
## #   age <dbl>, bh_cm <int>, bmi <dbl>, bsa_m2 <dbl>, lbm_kg <dbl>

C_0_0_lm2_v1 <- lm(log(C_0_0_order) ~ study_id + sex+age , data = order_0_df)
summary(C_0_0_lm2_v1) #sex and age depend on study_id, there is no point in looking at all the data

## 
## Call:
## lm(formula = log(C_0_0_order) ~ study_id + sex + age, data = order_0_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27793 -0.09038  0.04045  0.09746  0.15628 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)               -0.139503   0.141562  -0.985    0.341
## study_idWilkinson-1976-PO  0.052919   0.119010   0.445    0.663
## sexM                       0.025163   0.096885   0.260    0.799
## age                        0.001099   0.003004   0.366    0.720
## 
## Residual standard error: 0.1462 on 14 degrees of freedom
##   (128 пропущенных наблюдений удалены)
## Multiple R-squared:  0.04178,    Adjusted R-squared:  -0.1635 
## F-statistic: 0.2035 on 3 and 14 DF,  p-value: 0.8923

df_sa <- order_0_df %>% filter(!is.na(sex), !is.na(age), !is.na(C_0_0_order)) 
# Wilkinson and Thierauf data, 18 observations

df_sa <- df_sa %>%  mutate(bmi = wt_kg / (bh_cm/100)^2, 
         bsa_m2 = sqrt(bh_cm * wt_kg / 3600)) # bsa fits better in the model

C_0_0_lm2_v2 <- lm(log(C_0_0_order) ~ log(dose_g_per_kg) + bsa_m2 + sex, data = df_sa) 
summary(C_0_0_lm2_v2)

## 
## Call:
## lm(formula = log(C_0_0_order) ~ log(dose_g_per_kg) + bsa_m2 + 
##     sex, data = df_sa)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.24681 -0.08761  0.01368  0.09441  0.20278 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)  
## (Intercept)         -0.3044     0.3613  -0.842   0.4137  
## log(dose_g_per_kg)   0.9798     0.6431   1.524   0.1499  
## bsa_m2               0.4939     0.2732   1.808   0.0922 .
## sexM                -0.2210     0.1658  -1.333   0.2039  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.132 on 14 degrees of freedom
## Multiple R-squared:  0.2189, Adjusted R-squared:  0.05151 
## F-statistic: 1.308 on 3 and 14 DF,  p-value: 0.311

check_model(C_0_0_lm2_v2)

It looks better visually, but there are a lot of predictors. I calculated BMI and BSA and tried to include them in the model in various ways.

C_0_0_lm3_v1 <- lm(log(C_0_0_order) ~ bsa_m2 , data = order_0_df)

summary(C_0_0_lm3_v1)

## 
## Call:
## lm(formula = log(C_0_0_order) ~ bsa_m2, data = order_0_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25349 -0.10107  0.02324  0.09509  0.17215 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.4567     0.3092  -1.477    0.159
## bsa_m2        0.2087     0.1611   1.296    0.213
## 
## Residual standard error: 0.1329 on 16 degrees of freedom
##   (128 пропущенных наблюдений удалены)
## Multiple R-squared:  0.09499,    Adjusted R-squared:  0.03842 
## F-statistic: 1.679 on 1 and 16 DF,  p-value: 0.2134

check_model(C_0_0_lm3_v1)

1.6.5. Regression models for Vd

fit_Vd_1 <- lm(V_d_0_order ~ route + food + site, data = order_0_df_factor)

summary(fit_Vd_1)

## 
## Call:
## lm(formula = V_d_0_order ~ route + food + site, data = order_0_df_factor)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.187399 -0.045232  0.001478  0.044736  0.216100 
## 
## Coefficients: (1 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.61519    0.02098  29.320  < 2e-16 ***
## routePO      0.06870    0.02069   3.320  0.00129 ** 
## foodfed     -0.01162    0.01688  -0.688  0.49290    
## sitevenous        NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07512 on 93 degrees of freedom
##   (50 пропущенных наблюдений удалены)
## Multiple R-squared:  0.1337, Adjusted R-squared:  0.115 
## F-statistic: 7.174 on 2 and 93 DF,  p-value: 0.001266

cat("\n===============================================================\n")

## 
## ===============================================================

fit_Vd_3 <- lm(V_d_0_order ~ route, data = order_0_df_factor)

summary(fit_Vd_3)

## 
## Call:
## lm(formula = V_d_0_order ~ route, data = order_0_df_factor)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.195744 -0.041954  0.002783  0.046007  0.207755 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.60744    0.01766  34.402  < 2e-16 ***
## routePO      0.07317    0.01959   3.735 0.000322 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07491 on 94 degrees of freedom
##   (50 пропущенных наблюдений удалены)
## Multiple R-squared:  0.1292, Adjusted R-squared:   0.12 
## F-statistic: 13.95 on 1 and 94 DF,  p-value: 0.0003215

cat("\n===============================================================\n")

## 
## ===============================================================

fit_Vd_4 <- lm(V_d_0_order ~ route + site, data = order_0_df_factor)

summary(fit_Vd_4)

## 
## Call:
## lm(formula = V_d_0_order ~ route + site, data = order_0_df_factor)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.187399 -0.045232  0.001478  0.044736  0.216100 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.61519    0.02098  29.320  < 2e-16 ***
## routePO      0.06870    0.02069   3.320  0.00129 ** 
## sitevenous  -0.01162    0.01688  -0.688  0.49290    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07512 on 93 degrees of freedom
##   (50 пропущенных наблюдений удалены)
## Multiple R-squared:  0.1337, Adjusted R-squared:  0.115 
## F-statistic: 7.174 on 2 and 93 DF,  p-value: 0.001266

cat("\n===============================================================\n")

## 
## ===============================================================

fit_Vd_5 <- lm(V_d_0_order ~ route + food, data = order_0_df_factor)

summary(fit_Vd_5)

## 
## Call:
## lm(formula = V_d_0_order ~ route + food, data = order_0_df_factor)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.187399 -0.045232  0.001478  0.044736  0.216100 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.61519    0.02098  29.320  < 2e-16 ***
## routePO      0.06870    0.02069   3.320  0.00129 ** 
## foodfed     -0.01162    0.01688  -0.688  0.49290    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07512 on 93 degrees of freedom
##   (50 пропущенных наблюдений удалены)
## Multiple R-squared:  0.1337, Adjusted R-squared:  0.115 
## F-statistic: 7.174 on 2 and 93 DF,  p-value: 0.001266

cat("\n===============================================================\n")

## 
## ===============================================================

anova(fit_Vd_1, fit_Vd_3)

## Analysis of Variance Table
## 
## Model 1: V_d_0_order ~ route + food + site
## Model 2: V_d_0_order ~ route
##   Res.Df     RSS Df  Sum of Sq      F Pr(>F)
## 1     93 0.52484                            
## 2     94 0.52752 -1 -0.0026745 0.4739 0.4929

anova(fit_Vd_1, fit_Vd_4)

## Analysis of Variance Table
## 
## Model 1: V_d_0_order ~ route + food + site
## Model 2: V_d_0_order ~ route + site
##   Res.Df     RSS Df Sum of Sq F Pr(>F)
## 1     93 0.52484                      
## 2     93 0.52484  0         0

anova(fit_Vd_1, fit_Vd_5)

## Analysis of Variance Table
## 
## Model 1: V_d_0_order ~ route + food + site
## Model 2: V_d_0_order ~ route + food
##   Res.Df     RSS Df Sum of Sq F Pr(>F)
## 1     93 0.52484                      
## 2     93 0.52484  0         0

fit_Vd_3_sub <- lm(V_d_0_order ~ 1,  data = order_0_df_factor)

anova(fit_Vd_3, fit_Vd_3_sub)

## Analysis of Variance Table
## 
## Model 1: V_d_0_order ~ route
## Model 2: V_d_0_order ~ 1
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     94 0.52752                                  
## 2     95 0.60581 -1 -0.078296 13.952 0.0003215 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_model(fit_Vd_3)

fit_Vd_3_add <- lm(V_d_0_order ~ route + wt_kg * sex, data = order_0_df_factor)
fit_Vd_3_add2 <- lm(V_d_0_order ~ route + wt_kg + sex, data = order_0_df_factor)

anova(fit_Vd_3_add, fit_Vd_3_add2)

## Analysis of Variance Table
## 
## Model 1: V_d_0_order ~ route + wt_kg * sex
## Model 2: V_d_0_order ~ route + wt_kg + sex
##   Res.Df     RSS Df   Sum of Sq     F Pr(>F)
## 1     19 0.12279                            
## 2     20 0.12279 -1 -1.1327e-06 2e-04 0.9896

summary(fit_Vd_3_add2)

## 
## Call:
## lm(formula = V_d_0_order ~ route + wt_kg + sex, data = order_0_df_factor)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.120547 -0.051799 -0.007649  0.039997  0.182814 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.817074   0.116676   7.003 8.55e-07 ***
## routePO      0.016822   0.038723   0.434   0.6686    
## wt_kg       -0.003918   0.001551  -2.526   0.0201 *  
## sexM         0.110812   0.045764   2.421   0.0251 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07836 on 20 degrees of freedom
##   (122 пропущенных наблюдений удалены)
## Multiple R-squared:  0.2998, Adjusted R-squared:  0.1948 
## F-statistic: 2.855 on 3 and 20 DF,  p-value: 0.06298

fit_Vd <- lm(V_d_0_order ~ wt_kg + sex, data = order_0_df_factor)

anova(fit_Vd, fit_Vd_3_add2)

## Analysis of Variance Table
## 
## Model 1: V_d_0_order ~ wt_kg + sex
## Model 2: V_d_0_order ~ route + wt_kg + sex
##   Res.Df     RSS Df Sum of Sq      F Pr(>F)
## 1     21 0.12395                           
## 2     20 0.12279  1 0.0011586 0.1887 0.6686

summary(fit_Vd)

## 
## Call:
## lm(formula = V_d_0_order ~ wt_kg + sex, data = order_0_df_factor)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.115225 -0.054629 -0.004847  0.036033  0.187856 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.836182   0.105960   7.891 1.02e-07 ***
## wt_kg       -0.003953   0.001519  -2.603   0.0166 *  
## sexM         0.105957   0.043513   2.435   0.0239 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07683 on 21 degrees of freedom
##   (122 пропущенных наблюдений удалены)
## Multiple R-squared:  0.2932, Adjusted R-squared:  0.2259 
## F-statistic: 4.356 on 2 and 21 DF,  p-value: 0.02615

check_model(fit_Vd)

resid_panel(fit_Vd, plots = c("lev", "cookd"))

ggpredict(fit_Vd, terms = c("wt_kg", "sex"), ci_level = 0.95) %>% plot(show_data = TRUE, log_y = FALSE, jitter = .01) +
  labs(title = "Predicted values of Vd",
    subtitle = "Adjusted R-squared:  0.2259",
    color = "Sex") +
  xlab("Weight (kg)") +
  ylab("Vd (L / kg)")

plot_model(fit_Vd, digits = 3, show.values = TRUE,  show.intercept = TRUE, value.offset = 0.3,  axis.labels = c("Weight (kg)", "Sex", "Intercept"), title = "Regression coeffecients with 95% CI") +
  theme_nice()

1.6.6. Regression models for Beta

fit_Beta_1 <- lm(Beta ~ dose_g_per_kg * route, data = order_0_df_factor)

summary(fit_Beta_1)

## 
## Call:
## lm(formula = Beta ~ dose_g_per_kg * route, data = order_0_df_factor)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.138e-03 -2.014e-04 -1.855e-05  1.265e-04  8.673e-04 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            0.0003570  0.0002436   1.465    0.146    
## dose_g_per_kg          0.0035761  0.0005877   6.085 2.64e-08 ***
## routePO                0.0044200  0.0005433   8.135 1.88e-12 ***
## dose_g_per_kg:routePO -0.0073689  0.0009400  -7.840 7.76e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0003303 on 92 degrees of freedom
##   (50 пропущенных наблюдений удалены)
## Multiple R-squared:  0.5181, Adjusted R-squared:  0.5024 
## F-statistic: 32.97 on 3 and 92 DF,  p-value: 1.466e-14

fit_Beta_1_sub <- lm(Beta ~ 1,  data = order_0_df_factor)

anova(fit_Beta_1, fit_Beta_1_sub)

## Analysis of Variance Table
## 
## Model 1: Beta ~ dose_g_per_kg * route
## Model 2: Beta ~ 1
##   Res.Df        RSS Df   Sum of Sq      F    Pr(>F)    
## 1     92 1.0039e-05                                    
## 2     95 2.0830e-05 -3 -1.0791e-05 32.966 1.466e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_model(fit_Beta_1)

resid_panel(fit_Beta_1, plots = c("lev", "cookd"))

ggpredict(fit_Beta_1, terms = c("dose_g_per_kg", "route"), ci_level = 0) %>% plot(show_data = TRUE, log_y = FALSE, jitter = .01) + labs(subtitle = "Adjusted R-squared:  0.5671") + xlab("Dose (g per kg)")

plot_model(fit_Beta_1, show.values = TRUE,  show.intercept = TRUE, value.offset = 0.3,  axis.labels = c("Dose (g per kg)", "Route (PO)", "Dose (g per kg) x Route (PO)", "Intercept"), title = "Regression coeffecients with 95% CI") +
  theme_nice()

1.7. Interpretation of models. Effect of factors on PK-parameters

1.7.1. The effect of body mass and body surface area on PK-parameters

Body mass

ggpredict(fit_Vd, terms = c("wt_kg", "sex"), ci_level = 0.95) %>% plot(show_data = TRUE, log_y = FALSE, jitter = .01) +
  labs(title = "Predicted values of Vd",
    subtitle = "Adjusted R-squared:  0.2259",
    color = "Sex") +
  xlab("Weight (kg)") +
  ylab("Vd (L / kg)")

fit_Vd %>% broom::tidy(conf.int = TRUE) %>% filter(term == "wt_kg") %>% 
  mutate(across(c(where(is.numeric), -p.value), ~round(.x, 3))) %>% 
  mutate(p.value = gtsummary::style_pvalue(p.value),
         term = "Body mass, kg") %>% 
  select(-statistic) %>% 
  slice_tail(n = 1) %>% 
  rename(
    `Indicator` = term,
    `Coefficient (β)` = estimate,
    `Standard error` = std.error,
    `P-value` = p.value,
    `Lower bound of 95% CI` = conf.low,
    `Upper bound of 95% CI` = conf.high
  ) %>% 
  flextable() %>%
  add_header_row(., values = "Estimation of the body mass coefficient in the Vd model", colwidth = 6) %>% 
  flex_default_fun()

Estimation of the body mass coefficient in the Vd model
Indicator	Coefficient (β)	Standard error	P-value	Lower bound of 95% CI	Upper bound of 95% CI
Body mass, kg	-0.004	0.002	0.017	-0.007	-0.001

Body surface area

ggpredict(lm(V_d_0_order ~ bsa_m2 + sex, data = order_0_df),
          terms = c("bsa_m2", "sex"), ci_level = 0.95) %>% plot(show_data = TRUE, log_y = FALSE, jitter = .01) +
  labs(title = "Predicted values of Vd",
    subtitle = "Adjusted R-squared:  0.2433",
    color = "Sex") +
  xlab("Body Surface Area (m²)") +
  ylab("Vd (L / kg)")

lm(V_d_0_order ~ bsa_m2 + sex, data = order_0_df) %>% 
  broom::tidy(conf.int = TRUE) %>% filter(term == "bsa_m2") %>% 
  mutate(across(c(where(is.numeric), -p.value), ~round(.x, 3))) %>% 
  mutate(p.value = gtsummary::style_pvalue(p.value),
         term = "Body surface area, m²") %>% 
  select(-statistic) %>% 
  slice_tail(n = 1) %>% 
  rename(
    `Indicator` = term,
    `Coefficient (β)` = estimate,
    `Standard error` = std.error,
    `P-value` = p.value,
    `Lower bound of 95% CI` = conf.low,
    `Upper bound of 95% CI` = conf.high
  ) %>% 
  flextable() %>%
  add_header_row(., values = "Estimation of the body surface area coefficient in the Vd model", colwidth = 6) %>% 
  flex_default_fun()

Estimation of the body surface area coefficient in the Vd model
Indicator	Coefficient (β)	Standard error	P-value	Lower bound of 95% CI	Upper bound of 95% CI
Body surface area, m²	-0.311	0.131	0.031	-0.589	-0.032

Lean body mass

ggpredict(Cmax_lm2_v3,
          terms = c("lbm_kg", "age"), ci_level = 0.95) %>% plot(show_data = TRUE, log_y = FALSE, jitter = .01) +
  labs(title = "Predicted values of Cmax",
    subtitle = "Adjusted R-squared:  0.3328",
    color = "Age") +
  xlab("Lean Body Mass (kg)") +
  ylab("Cmax (g/L)")

Cmax_lm2_v3 %>% 
  broom::tidy(conf.int = TRUE) %>%
  # filter(term == "lbm_kg") %>% 
  mutate(across(c(where(is.numeric), -p.value), ~round(.x, 4))) %>% 
  mutate(p.value = gtsummary::style_pvalue(p.value),
         term = c("Intersept", "Age, years", "Lean body mass, kg")) %>% 
  select(-statistic) %>% 
  rename(
    `Indicator` = term,
    `Coefficient (β)` = estimate,
    `Standard error` = std.error,
    `P-value` = p.value,
    `Lower bound of 95% CI` = conf.low,
    `Upper bound of 95% CI` = conf.high
  ) %>% 
  flextable() %>%
  add_header_row(., values = "Estimation of the lean body mass coefficient in the Cmax model", colwidth = 6) %>% 
  flex_default_fun()

Estimation of the lean body mass coefficient in the Cmax model
Indicator	Coefficient (β)	Standard error	P-value	Lower bound of 95% CI	Upper bound of 95% CI
Intersept	1.8629	0.3048	<0.001	1.2498	2.4760
Age, years	0.0053	0.0019	0.008	0.0015	0.0091
Lean body mass, kg	-0.0086	0.0042	0.045	-0.0170	-0.0002

With an increase in body weight by 1 kg, Vd decreases by 0.004 (95% CI [0.001; 0.007]) L/kg (all other things being equal), and with an increase in body surface area by 1 \(m^2\), Vd decreases by 0.311 (95% CI [0.032; 0.589]) L/kg (all other things being equal).

When lean body mass increases by 1 kg, Cmax decreases by 0.0086 (95% CI [0.0002; 0.017]) g/L, all other things being equal.

1.7.2. The effect of sex on PK-parameters

Sex does not significantly affect AUC (Male 95%CI[-0.3108; 0.1754]) or Cmax (Male 95%CI[-0.1505; 0.326]).

In male subjects, the volume of distribution is 0.106 L/kg higher than in females, all other things being equal (95% CI [-0.0155; 0.1964]).

fit_Vd %>% 
  broom::tidy(conf.int = TRUE) %>% filter(term == "sexM") %>% 
  mutate(across(c(where(is.numeric), -p.value), ~round(.x, 4))) %>% 
  mutate(p.value = gtsummary::style_pvalue(p.value),
         term = "Sex") %>% 
  select(-statistic) %>% 
  slice_tail(n = 1) %>% 
  rename(
    `Indicator` = term,
    `Coefficient (β)` = estimate,
    `Standard error` = std.error,
    `P-value` = p.value,
    `Lower bound of 95% CI` = conf.low,
    `Upper bound of 95% CI` = conf.high
  ) %>% 
  flextable() %>%
  add_header_row(., values = " Estimation of the sex coefficient in the Vd model", colwidth = 6) %>% 
  flex_default_fun()

Estimation of the sex coefficient in the Vd model
Indicator	Coefficient (β)	Standard error	P-value	Lower bound of 95% CI	Upper bound of 95% CI
Sex	0.106	0.0435	0.024	0.0155	0.1964

1.7.3. The effect of food on PK-parameters

It’s necessary to replace the site to food from the previous model of AUC (fit_AUC_mod_int), because you can’t add them together.

fit_AUC_mod_with_food <- lm(log(AUC) ~ dose_g_per_kg*route + food, order_0_df)

summary(fit_AUC_mod_with_food, confint = TRUE, digits = 3)

## 
## Call:
## lm(formula = log(AUC) ~ dose_g_per_kg * route + food, data = order_0_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47726 -0.11435  0.01464  0.09105  0.55281 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             4.5202     0.1594  28.365  < 2e-16 ***
## dose_g_per_kg           0.6720     0.3498   1.921 0.057871 .  
## routePO                -1.0788     0.2868  -3.762 0.000298 ***
## foodfed                -0.3576     0.0459  -7.791 1.04e-11 ***
## dose_g_per_kg:routePO   2.1069     0.4977   4.233 5.50e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1744 on 91 degrees of freedom
##   (50 пропущенных наблюдений удалены)
## Multiple R-squared:  0.8079, Adjusted R-squared:  0.7995 
## F-statistic: 95.68 on 4 and 91 DF,  p-value: < 2.2e-16

plot_auc <- plot_model(fit_AUC_mod_with_food, show.values = TRUE,  show.intercept = TRUE, value.offset = 0.3, axis.labels = c("Dose (g per kg) x Route (PO)", "Food status (fed)", "Route (PO)", "Dose (g per kg)","Intercept")) +
  geom_hline(yintercept = 0, color = "gray3", size = 0.5) +
  labs(title = "Regression coefficients for AUC with 95% CI", 
       subtitle = "log(AUC) ~ dose_g_per_kg*route + food") +
  theme_nice()+
  theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

plot_auc

d <- plot_auc$data

d$ci_label <- sprintf(
  "[%.2f; %.2f]",
  d$conf.low,
  d$conf.high
)

Compared with the fasted state, “fed” condition was associated with log(AUC) decreased by 0.36, AUC decreased by 30.2% (all other things being equal).Percent derived as (exp(β)−1)×100.

plot_cmax <- plot_model(Cmax_lm2_v2, show.values = TRUE,  show.intercept = TRUE, value.offset = 0.3, digits = 3) +
  geom_hline(yintercept = 0, color = "gray3", size = 0.5) +
  labs(title = "Regression coefficients for Cmax with 95% CI", 
       subtitle = "log(Cmax) ~ dose_g_per_kg + route + food + sex + age") +
  theme_nice()+
  theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

plot_cmax

Compared with the fasted state, “fed” condition was associated with log(Cmax) decreased by 0.286, Cmax decreased by 24.9% (all other things being equal).

plot_beta <- plot_model(fit_Beta_1, show.values = TRUE,  show.intercept = TRUE, value.offset = 0.5, digits = 3,   axis.labels = c("Dose (g per kg) x Route (PO)", "food (fed)","Route (PO)", "Dose (g per kg)", "Intercept")) +
    geom_hline(yintercept = 0, color = "gray3", size = 0.5) +
labs(title = "Regression coefficients \n for Beta with 95% CI", 
       subtitle = "Beta ~ dose_g_per_kg * route + food") +
  scale_y_continuous(limits = c(-0.015, 0.015)) +
  theme_nice()+
  theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

plot_beta

d <- plot_beta$data

d$ci_label <- sprintf(
  "[%.4f; %.4f]",
  d$conf.low,
  d$conf.high)

Compared with the fasted state, “fed” condition was associated with Beta increased by 0.0004 g/L/min (all other things being equal).

2. Allometric scaling using rats’ PK parameters

2.1. Data unification

Kozawa_2007 <- Kozawa_2007 %>% 
  mutate(time_min = round(time_h * 60, 0)) %>% 
  mutate(route = route %>% str_remove_all(" bolus")) %>%
  select(-time_h) %>% 
  mutate(number = dose_g_per_kg/0.5)

Eriksson_1977 <- Eriksson_1977 %>% 
  mutate(time_min = round(time_h * 60, 0)) %>% 
  select(-time_h) %>% 
  mutate(number = case_when(
      dose_g_per_kg == 0.75 ~ 1,
      dose_g_per_kg == 1.50 ~ 2,
      dose_g_per_kg == 3.00 ~ 3,
    ))

2.2. Concentration-time plots

ggplot(Kozawa_2007, aes(x = time_min, y = conc_g_per_L)) +
  geom_line(aes(group = number)) +
  scale_y_continuous(breaks = seq(0, 3.6, 0.2)) +
  scale_x_continuous(breaks = seq(0, 500, 25)) +
  labs(title = "Concentration-time plot (rats, IV route)",
       subtitle = "Averaged data",
       caption = "Source: Kozawa-2007",
       x = "Time, min",
       y = "Concentration, g/L")

ggplot(Eriksson_1977, aes(x = time_min, y = conc_g_per_L)) +
  geom_line(aes(group = number)) +
  scale_y_continuous(breaks = seq(0, 1.8, 0.1)) +
  scale_x_continuous(breaks = seq(0, 300, 30)) +
  labs(title = "Concentration-time plot (rats, PO route)",
       subtitle = "Averaged data",
       caption = "Source: Eriksson-1977",
       x = "Time, min",
       y = "Concentration, g/L")

2.3. PK-parameteres calculation

Kozawa_2007_pk <- Kozawa_2007 %>% 
  arrange(time_min) %>% 
  nest(.by = number, .key = "average_data") %>% 
  mutate(pharma = purrr::map(average_data,
                            ~pharma_fun(.$conc_g_per_L,
                                        .$time_min,
                                        .$dose_g_per_kg))) %>% 
  unnest_wider(col = pharma) %>% 
  select(-average_data)

Kozawa_2007_pk <- Kozawa_2007 %>% 
  select(-c(time_min, conc_g_per_L)) %>% 
  group_by(number) %>% 
  slice(1) %>% 
  ungroup() %>% 
  right_join(Kozawa_2007_pk) %>% 
  select(-number)

## Joining with `by = join_by(number)`

Eriksson_1977_pk <- Eriksson_1977 %>% 
  arrange(time_min) %>% 
  nest(.by = number, .key = "average_data") %>% 
  mutate(pharma = purrr::map(average_data,
                            ~pharma_fun(.$conc_g_per_L,
                                        .$time_min,
                                        .$dose_g_per_kg))) %>% 
  unnest_wider(col = pharma) %>% 
  select(-average_data)

# for two concentrations there are not enough points to calculate beta

Eriksson_1977_pk <- Eriksson_1977 %>% 
  select(-c(time_min, conc_g_per_L)) %>% 
  group_by(number) %>% 
  slice(1) %>% 
  ungroup() %>% 
  right_join(Eriksson_1977_pk) %>% 
  select(-number)

## Joining with `by = join_by(number)`

rat_data_pk <- rbind(Kozawa_2007_pk, Eriksson_1977_pk)

2.4. Preaparing data and function

# Dataset prepararion
rat_data_pk <-  rat_data_pk %>% 
  mutate(V_max_0_order = Beta * 60 * V_d_0_order,
         BW_g = BW_g / 1000) %>%
  rename(wt_kg = BW_g) %>% 
  relocate(V_max_0_order, .after = V_d_0_order)

allom_scale <- function(V_rat, # L, abs
                        rat_weight, # kg
                        human_weight = 1, # kg
                        n = 1)
{
  V_human = V_rat * (human_weight / rat_weight) ^ n
  return(V_human)
}

# Human dataset preparation
allometric <- order_0_df %>% 
  filter(!is.na(wt_kg)) %>% 
  select(study_id, subject_id, dose_g_per_kg, route, 
         Beta, V_d_0_order, wt_kg) %>% 
  mutate(V_max = Beta * V_d_0_order * 60) %>% 
  mutate(Vd_obs = V_d_0_order * wt_kg,
         Vmax_obs = Beta * V_d_0_order * wt_kg * 60,
         Beta_obs = Beta * 60) %>% 
  relocate(c(Vd_obs, Vmax_obs, Beta_obs), .after = V_max)

# Only Wilkinson IV (Human dose 0,5-0,62 mg/kg)
# allometric_IV_Wilk <- allometric %>% 
#   filter(route == "IV",
#          study_id == "Wilkinson-1976-IV") %>% 
#   crossing(rat_data_pk %>% filter(dose_g_per_kg == 0.5) %>% 
#              mutate(Vd_abs_rat = (V_d_0_order * wt_kg)) %>% 
#              rename_with(~ paste0("rat_", .x))) %>% 
#   mutate(Vd_classic = allom_scale(rat_Vd_abs_rat, rat_wt_kg, wt_kg),
#          Vd_matsumoto = allom_scale(rat_Vd_abs_rat, rat_wt_kg, wt_kg, 0.932),
#          Beta_classic = allom_scale(rat_Beta, rat_wt_kg, wt_kg, 0.75),
#          Beta_matsumoto = allom_scale(rat_Beta, rat_wt_kg, wt_kg, 0.708),
#          Vmax_classic = allom_scale(rat_V_max_0_order, rat_wt_kg, wt_kg, 0.75),
#          Vmax_matsumoto = allom_scale(rat_V_max_0_order, rat_wt_kg, wt_kg, 0.708))

2.5. Allometric scaling

2.5.1. IV allometric scaling (main df)

# Allometric data for IV route
allometric_IV <- allometric %>% 
  filter(route == "IV") %>% 
  crossing(rat_data_pk %>% filter(dose_g_per_kg == 0.5) %>% 
             mutate(Vd_abs_rat = (V_d_0_order * wt_kg)) %>% 
             rename_with(~ paste0("rat_", .x))) %>% 
  mutate(Vd_classic = allom_scale(rat_Vd_abs_rat, rat_wt_kg, wt_kg),
         Vd_matsumoto = allom_scale(rat_Vd_abs_rat, rat_wt_kg, wt_kg, 0.932),
         Beta_classic = allom_scale(rat_Beta * 60, rat_wt_kg, wt_kg, 0.75 - 1),
         Beta_matsumoto = allom_scale(rat_Beta * 60, rat_wt_kg, wt_kg, 0.708 - 0.932),
         Vmax_classic = allom_scale(rat_V_max_0_order, rat_wt_kg, wt_kg, 0.75),
         Vmax_matsumoto = allom_scale(rat_V_max_0_order, rat_wt_kg, wt_kg, 0.708)) %>% 
  mutate(Vd_metr_cl = (Vd_classic - Vd_obs) / Vd_obs * 100,
         Vd_metr_mats = (Vd_matsumoto - Vd_obs) / Vd_obs * 100,
         Beta_metr_cl = (Beta_classic - Beta_obs) / Beta_obs * 100,
         Beta_metr_mats = (Beta_matsumoto - Beta_obs) / Beta_obs * 100,
         Vmax_metr_cl = (Vmax_classic - Vmax_obs) / Vmax_obs * 100,
         Vmax_metr_mats = (Vmax_matsumoto - Vmax_obs) / Vmax_obs * 100
  ) %>% 
  mutate(across(contains("metr"), ~round(.x, 3)))

# Observed and mean dose
main_df_IV <- order_0_df %>% 
  filter(route == "IV" & !is.na(wt_kg)) %>% 
  mutate(Vd_obs = V_d_0_order * wt_kg,
         Vmax_obs = Beta * V_d_0_order * wt_kg * 60,
         Beta_obs = Beta * 60)

cat("Human mean dose, IV route, g/kg:", mean(main_df_IV$dose_g_per_kg))

## Human mean dose, IV route, g/kg: 0.5708929

allometric_IV %>% 
  select(contains("Vd_metr")) %>% 
  rename(`Allometric classic` = Vd_metr_cl,
         `Matsumoto` = Vd_metr_mats) %>% 
  pivot_longer(everything(), names_to = "Type", values_to = "Percentage") %>% 
  filter(!is.na(Percentage)) %>% 
ggplot(aes(x = Type, y = Percentage)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  labs(
    x = "Used allometric scale exponent",
    y = "Error, %", 
    title = "Error of Vd allometric scaling using different allometric\nscale exponents, IV route of administration",
    subtitle = "Rat dose = 0.5 g/kg"
  ) +
  theme_bw()

# Vd
Vd_vis_IV <- bind_rows(
  main_df_IV %>% transmute(value = Vd_obs, model = "Observed"),
  allometric_IV %>% transmute(value = Vd_classic, model = "Allometric classic"),
  allometric_IV %>% transmute(value = Vd_matsumoto, model = "Matsumoto")
)

ggplot(Vd_vis_IV, aes(x = model, y = value)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  theme_bw()+
  labs(
    x = "Used allometric scale exponent",
    y = "Vd, L", 
    title = "Vd, allometric scaling using different allometric\nscale exponents, IV route of administration",
    subtitle = "Mean human dose = 0.57 g/kg, rat dose = 0.5 g/kg"
  )

## Warning: Removed 50 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

The usage of classical exponent (=1) for Vd allometric scaling shows mean error (%) from +20 to +50%. Matsumoto’s exponent (=0,932) shows mean error (%) from 0% to -15%.

allometric_IV %>% 
  select(contains("Beta_metr")) %>% 
  rename(`Allometric classic` = Beta_metr_cl,
         `Matsumoto` = Beta_metr_mats) %>% 
  pivot_longer(everything(), names_to = "Type", values_to = "Percentage") %>% 
  filter(!is.na(Percentage)) %>% 
ggplot(aes(x = Type, y = Percentage)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  labs(
    x = "Used allometric scale exponent",
    y = "Error, %", 
    title = "Error of Beta allometric scaling using different allometric\nscale exponents, IV route of administration",
    subtitle = "Rat dose = 0.5 g/kg"
  ) +
  theme_bw()

# Beta
Beta_vis_IV <- bind_rows(
  main_df_IV %>% transmute(value = Beta_obs, model = "Observed"),
  allometric_IV %>% transmute(value = Beta_classic, model = "Allometric classic"),
  allometric_IV %>% transmute(value = Beta_matsumoto, model = "Matsumoto")
)

ggplot(Beta_vis_IV, aes(x = model, y = value)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  theme_bw()+
  labs(
    x = "Used allometric scale exponent",
    y = "Beta, g/L/hour", 
    title = "Beta, allometric scaling using different allometric\nscale exponents, IV route of administration",
    subtitle = "Mean human dose = 0.57 g/kg, rat dose = 0.5 g/kg"
  )

## Warning: Removed 50 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

The usage of classical exponent (= 0,75 - 1) for Beta allometric scaling shows mean error (%) from -35% to -50%. Matsumoto’s exponent (= 0,708 - 0,932) shows mean error (%) from -25% to -43%.

allometric_IV %>% 
  select(contains("Vmax_metr")) %>% 
  rename(`Allometric classic` = Vmax_metr_cl,
         `Matsumoto` = Vmax_metr_mats) %>% 
  pivot_longer(everything(), names_to = "Type", values_to = "Percentage") %>% 
  filter(!is.na(Percentage)) %>% 
ggplot(aes(x = Type, y = Percentage)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  labs(
    x = "Used allometric scale exponent",
    y = "Error, %", 
    title = "Error of Vmax allometric scaling using different allometric\nscale exponents, IV route of administration",
    subtitle = "Rat dose = 0.5 g/kg"
  ) +
  theme_bw()

# Vmax
Vmax_vis_IV <- bind_rows(
  main_df_IV %>% transmute(value = Vmax_obs, model = "Observed"),
  allometric_IV %>% transmute(value = Vmax_classic, model = "Allometric classic"),
  allometric_IV %>% transmute(value = Vmax_matsumoto, model = "Matsumoto")
)

ggplot(Vmax_vis_IV, aes(x = model, y = value)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  theme_bw()+
  labs(
    x = "Used allometric scale exponent",
    y = "Vmax, g/hour", 
    title = "Vmax, allometric scaling using different allometric\nscale exponents, IV route of administration",
    subtitle = "Mean human dose = 0.57 g/kg, rat dose = 0.5 g/kg"
  )

## Warning: Removed 50 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

The usage of classical exponent (= 0,75) for Vmax allometric scaling shows mean error (%) from +87 to +150%. Matsumoto’s exponent (= 0,708) shows mean error (%) from +50% to +100%.

2.5.2. PO allometric scaling

2.5.2.1. PO allometric scaling (750 mg per kg)

# Allometric scaling
allometric_PO <- allometric %>% 
  filter(route == "PO") %>% 
  crossing(rat_data_pk %>% filter(dose_g_per_kg == 0.75) %>% 
             mutate(Vd_abs_rat = (V_d_0_order * wt_kg)) %>% 
             rename_with(~ paste0("rat_", .x))) %>% 
  mutate(Vd_classic = allom_scale(rat_Vd_abs_rat, rat_wt_kg, wt_kg),
         Vd_matsumoto = allom_scale(rat_Vd_abs_rat, rat_wt_kg, wt_kg, 0.932)) %>% 
  mutate(Vd_metr_cl = (Vd_classic - Vd_obs) / Vd_obs * 100,
         Vd_metr_mats = (Vd_matsumoto - Vd_obs) / Vd_obs * 100,
  ) %>% 
  mutate(across(contains("metr"), ~round(.x, 3)))

# Human dataset preparation
main_df_PO <- order_0_df %>% 
  filter(route == "PO" & !is.na(wt_kg)) %>% 
  mutate(Vd_obs = V_d_0_order * wt_kg,
         Vmax_obs = Beta * V_d_0_order * wt_kg * 60,
         Beta_obs = Beta * 60)

cat("Human mean dose, PO route, g/kg):", mean(main_df_PO$dose_g_per_kg))

## Human mean dose, PO route, g/kg): 0.5801667

allometric_PO %>% 
  select(contains("Vd_metr")) %>% 
  rename(`Allometric classic` = Vd_metr_cl,
         `Matsumoto` = Vd_metr_mats) %>% 
  pivot_longer(everything(), names_to = "Type", values_to = "Percentage") %>% 
  filter(!is.na(Percentage)) %>% 
ggplot(aes(x = Type, y = Percentage)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  labs(
    x = "Used allometric scale exponent",
    y = "Error, %", 
    title = "Error of Vd allometric scaling using different allometric\nscale exponents, PO route of administration",
    subtitle = "Rat dose = 0.75 g/kg"
  ) +
  theme_bw()

Mean error (%) for Vd allometric scaling is more than +100% up to +300% using classic allometric exponent. We cannot use another PK parameters due to no captured elimination phase in the rat data for this dose.

# Vd
Vd_vis_PO <- bind_rows(
  main_df_PO %>% transmute(value = Vd_obs, model = "Observed"),
  allometric_PO %>% transmute(value = Vd_classic, model = "Allometric classic"),
  allometric_PO %>% transmute(value = Vd_matsumoto, model = "Matsumoto")
)

ggplot(Vd_vis_PO, aes(x = model, y = value)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  theme_bw()+
  labs(
    x = "Used allometric scale exponent",
    y = "Vd, L", 
    title = "Vd, allometric scaling using different allometric\nscale exponents, PO route of administration",
    subtitle = "Mean human dose = 0.58 g/kg, rat dose = 0.75 g/kg"
  )

2.5.2.2. PO allometric scaling (3000 mg per kg)

# Scaling
allometric_PO_3 <- allometric %>% 
  filter(route == "PO") %>% 
  crossing(rat_data_pk %>% filter(dose_g_per_kg == 3) %>% 
             mutate(Vd_abs_rat = (V_d_0_order * wt_kg)) %>% 
             rename_with(~ paste0("rat_", .x))) %>%
  mutate(Vd_classic = allom_scale(rat_Vd_abs_rat, rat_wt_kg, wt_kg),
         Vd_matsumoto = allom_scale(rat_Vd_abs_rat, rat_wt_kg, wt_kg, 0.932),
         Beta_classic = allom_scale(rat_Beta * 60, rat_wt_kg, wt_kg, 0.75 - 1),
         Beta_matsumoto = allom_scale(rat_Beta * 60, rat_wt_kg, wt_kg, 0.708 - 0.932),
         Vmax_classic = allom_scale(rat_V_max_0_order, rat_wt_kg, wt_kg, 0.75),
         Vmax_matsumoto = allom_scale(rat_V_max_0_order, rat_wt_kg, wt_kg, 0.708)) %>% 
  mutate(Vd_metr_cl = (Vd_classic - Vd_obs) / Vd_obs * 100,
         Vd_metr_mats = (Vd_matsumoto - Vd_obs) / Vd_obs * 100,
         Beta_metr_cl = (Beta_classic - Beta_obs) / Beta_obs * 100,
         Beta_metr_mats = (Beta_matsumoto - Beta_obs) / Beta_obs * 100,
         Vmax_metr_cl = (Vmax_classic - Vmax_obs) / Vmax_obs * 100,
         Vmax_metr_mats = (Vmax_matsumoto - Vmax_obs) / Vmax_obs * 100
  ) %>% 
  mutate(across(contains("metr"), ~round(.x, 3)))

allometric_PO_3 %>% 
  select(contains("Vd_metr")) %>% 
  rename(`Allometric classic` = Vd_metr_cl,
         `Matsumoto` = Vd_metr_mats) %>% 
  pivot_longer(everything(), names_to = "Type", values_to = "Percentage") %>% 
  filter(!is.na(Percentage)) %>% 
ggplot(aes(x = Type, y = Percentage)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  labs(
    x = "Used allometric scale exponent",
    y = "Error, %", 
    title = "Error of Vd allometric scaling using different allometric\nscale exponents, PO route of administration",
    subtitle = "Rat dose = 3.00 g/kg"
  ) +
  theme_bw()

# Vd
Vd_vis_PO_3 <- bind_rows(
  main_df_PO %>% transmute(value = Vd_obs, model = "Observed"),
  allometric_PO_3 %>% transmute(value = Vd_classic, model = "Allometric classic"),
  allometric_PO_3 %>% transmute(value = Vd_matsumoto, model = "Matsumoto")
)

ggplot(Vd_vis_PO_3, aes(x = model, y = value)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  theme_bw()+
  labs(
    x = "Used allometric scale exponent",
    y = "Vd, L", 
    title = "Vd, allometric scaling using different allometric\nscale exponents, PO route of administration",
    subtitle = "Mean human dose = 0.58 g/kg, rat dose = 3.00 g/kg"
  )

allometric_PO_3 %>% 
  select(contains("Beta_metr")) %>% 
  rename(`Allometric classic` = Beta_metr_cl,
         `Matsumoto` = Beta_metr_mats) %>% 
  pivot_longer(everything(), names_to = "Type", values_to = "Percentage") %>% 
  filter(!is.na(Percentage)) %>% 
ggplot(aes(x = Type, y = Percentage)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  labs(
    x = "Used allometric scale exponent",
    y = "Error, %", 
    title = "Error of Beta allometric scaling using different allometric\nscale exponents, PO route of administration",
    subtitle = "Rat dose = 3.00 g/kg"
  ) +
  theme_bw()

# Beta
Beta_vis_PO_3 <- bind_rows(
  main_df_PO %>% transmute(value = Beta_obs, model = "Observed"),
  allometric_PO_3 %>% transmute(value = Beta_classic, model = "Allometric classic"),
  allometric_PO_3 %>% transmute(value = Beta_matsumoto, model = "Matsumoto")
)

ggplot(Beta_vis_PO_3, aes(x = model, y = value)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  theme_bw()+
  labs(
    x = "Used allometric scale exponent",
    y = "Beta, g/L/hour", 
    title = "Beta, allometric scaling using different allometric\nscale exponents, PO route of administration",
    subtitle = "Mean human dose = 0.58 g/kg, rat dose = 3.00 g/kg"
  )

allometric_PO_3 %>% 
  select(contains("Vmax_metr")) %>% 
  rename(`Allometric classic` = Vmax_metr_cl,
         `Matsumoto` = Vmax_metr_mats) %>% 
  pivot_longer(everything(), names_to = "Type", values_to = "Percentage") %>% 
  filter(!is.na(Percentage)) %>% 
ggplot(aes(x = Type, y = Percentage)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  labs(
    x = "Used allometric scale exponent",
    y = "Error, %", 
    title = "Error of Vmax allometric scaling using different allometric\nscale exponents, PO route of administration",
    subtitle = "Rat dose = 3.00 g/kg"
  ) +
  theme_bw()

# Vmax
Vmax_vis_PO_3 <- bind_rows(
  main_df_PO %>% transmute(value = Vmax_obs, model = "Observed"),
  allometric_PO_3 %>% transmute(value = Vmax_classic, model = "Allometric classic"),
  allometric_PO_3 %>% transmute(value = Vmax_matsumoto, model = "Matsumoto")
)

ggplot(Vmax_vis_PO_3, aes(x = model, y = value)) +
  geom_boxplot(width = 0.1, staplewidth = 0.5) +
  theme_bw()+
  labs(
    x = "Used allometric scale exponent",
    y = "Vmax, g/hour", 
    title = "Vmax, allometric scaling using different allometric\nscale exponents, PO route of administration",
    subtitle = "Mean human dose = 0.58 g/kg, rat dose = 3.00 g/kg"
  )

Predictive mean performance (error, %) for PO route with rat dose 3.0 g/kg:

• Vd: Allometric classic +187.82, Matsumoto +98.65;

• Beta: Allometric classic -72.20, Matsumoto -67.97;

• Vmax: Allometric classic +144.25, Matsumoto +94.28.